[orx-math] Add matrix and sparse matrix implementations with Cholesky and QR decompositions

orx-math/src/commonMain/kotlin/matrix/Matrix.kt

@@ -0,0 +1,177 @@
+package org.openrndr.extra.math.matrix
+
+/**
+ * Represents a two-dimensional matrix with support for basic operations such as indexing,
+ * copying, and mathematical computations.
+ *
+ * @property rows The number of rows in the matrix.
+ * @property cols The number of columns in the matrix.
+ */
+class Matrix(val rows: Int, val cols: Int) {
+    val data = Array(rows) { DoubleArray(cols) }
+
+    operator fun get(i: Int, j: Int) = data[i][j]
+    operator fun set(i: Int, j: Int, value: Double) {
+        data[i][j] = value
+    }
+
+    fun copy(): Matrix {
+        val result = Matrix(rows, cols)
+        for (i in 0 until rows) {
+            for (j in 0 until cols) {
+                result[i, j] = this[i, j]
+            }
+        }
+        return result
+    }
+
+    /**
+     * Returns a new matrix that is the transpose of this matrix.
+     *
+     * The transpose of a matrix is obtained by flipping the matrix over its diagonal,
+     * effectively switching the row and column indices of each element.
+     *
+     * @return A new matrix representing the transpose of the current matrix.
+     */
+    fun transposed() : Matrix {
+        val result = Matrix(cols, rows)
+        for (i in 0 until rows) {
+            for (j in 0 until cols) {
+                result[j, i] = this[i, j]
+            }
+        }
+        return result
+    }
+
+    /**
+     * Checks if the given matrix is symmetric within a specified tolerance.
+     *
+     * A matrix is considered symmetric if it is square (same number of rows and columns)
+     * and satisfies the condition matrix[i, j] == matrix[j, i] for all i and j, within
+     * the given tolerance.
+     *
+     * @param matrix The matrix to be checked for symmetry.
+     * @param tolerance The maximum allowable difference between matrix[i, j] and matrix[j, i]
+     *                  to still consider the matrix symmetric. Default is 1e-10.
+     * @return `true` if the matrix is symmetric within the specified tolerance, `false` otherwise.
+     */
+    fun isSymmetric(matrix: Matrix, tolerance: Double = 1e-10): Boolean {
+        if (matrix.rows != matrix.cols) return false
+
+        for (i in 0 until matrix.rows) {
+            for (j in 0 until matrix.cols) {
+                if (kotlin.math.abs(matrix[i, j] - matrix[j, i]) > tolerance) {
+                    return false
+                }
+            }
+        }
+        return true
+    }
+
+    companion object {
+        /**
+         * Generates an identity matrix of the specified size.
+         *
+         * An identity matrix is a square matrix with ones on the diagonal
+         * and zeros elsewhere.
+         *
+         * @param n The size of the identity matrix (number of rows and columns).
+         *          Must be a positive integer.
+         * @return A square matrix of size n x n, where all diagonal elements are 1
+         *         and all non-diagonal elements are 0.
+         */
+        fun identity(n: Int): Matrix {
+            val result = Matrix(n, n)
+            for (i in 0 until n) {
+                result[i, i] = 1.0
+            }
+            return result
+        }
+
+        /**
+         * Creates a new matrix with the specified number of rows and columns, filled with zeros.
+         *
+         * @param rows The number of rows in the matrix.
+         * @param cols The number of columns in the matrix.
+         * @return A matrix with dimensions [rows] x [cols], initialized to zero.
+         */
+        fun zeros(rows: Int, cols: Int): Matrix {
+            return Matrix(rows, cols)
+        }
+    }
+
+
+    /**
+     * Multiplies two matrices and returns the resulting matrix.
+     *
+     * @param a The first matrix to be multiplied.
+     * @param b The second matrix to be multiplied.
+     * @return The resulting matrix after multiplying A and B.
+     * @throws IllegalArgumentException If the number of columns in A does not match the number of rows in B.
+     */
+    fun multiply(a: Matrix, b: Matrix): Matrix {
+        if (a.cols != b.rows) {
+            throw IllegalArgumentException("Matrix dimensions don't match for multiplication")
+        }
+
+        val result = Matrix.zeros(a.rows, b.cols)
+        for (i in 0 until a.rows) {
+            for (j in 0 until b.cols) {
+                for (k in 0 until a.cols) {
+                    result[i, j] += a[i, k] * b[k, j]
+                }
+            }
+        }
+        return result
+    }
+
+    /**
+     * Multiplies this matrix with another matrix and returns the resulting matrix.
+     *
+     * @param other The matrix to multiply with this matrix.
+     * @return The resulting matrix after multiplication.
+     * @throws IllegalArgumentException If the number of columns in this matrix does not match the number of rows in the other matrix.
+     */
+    operator fun times(other: Matrix): Matrix {
+        return multiply(this, other)
+    }
+}
+
+fun Matrix.columnMean(): Matrix {
+    val means = DoubleArray(cols)
+
+    for (j in 0 until rows) {
+        for (i in 0 until cols) {
+            means[i] += this[j, i]
+        }
+
+    }
+    for (i in 0 until cols) {
+        means[i] /= rows.toDouble()
+    }
+    val result = Matrix.zeros(1, cols)
+    for (i in 0 until cols) {
+        result[0, i] = means[i]
+    }
+    return result
+}
+
+operator fun Matrix.minus(other: Matrix): Matrix {
+    val result = Matrix.zeros(rows, cols)
+    if (cols == other.cols && other.rows == 1) {
+        for (j in 0 until rows) {
+            for (i in 0 until cols) {
+                result[j, i] = this[j, i] - other[0, i]
+            }
+        }
+    } else if (cols == other.cols && rows == other.rows) {
+        for (j in 0 until rows) {
+            for (i in 0 until cols) {
+                result[j, i] = this[j, i] - other[j, i]
+            }
+        }
+    } else {
+        error("Cannot subtract matrices of different dimensions")
+    }
+    return result
+}

orx-math/src/commonMain/kotlin/matrix/MatrixCholesky.kt

@@ -0,0 +1,132 @@
+package org.openrndr.extra.math.matrix
+
+import kotlin.math.sqrt
+
+/**
+ * Performs the Cholesky decomposition on a given square, positive-definite matrix.
+ * The Cholesky decomposition expresses the matrix as L * L^T, where L is a lower triangular matrix.
+ *
+ * @param matrix The input matrix to decompose. Must be square (rows == cols) and positive-definite.
+ * @return A lower triangular matrix resulting from the Cholesky decomposition of the input matrix.
+ * @throws IllegalArgumentException if the input matrix is not square or not positive-definite.
+ */
+fun choleskyDecomposition(matrix: Matrix): Matrix {
+    if (matrix.rows != matrix.cols) {
+        throw IllegalArgumentException("Matrix must be square")
+    }
+
+    val n = matrix.rows
+    val L = Matrix.zeros(n, n)
+
+    for (i in 0 until n) {
+        for (j in 0..i) {
+            when {
+                i == j -> {
+                    // Diagonal elements
+                    var sum = 0.0
+                    for (k in 0 until j) {
+                        sum += L[j, k] * L[j, k]
+                    }
+                    val diagonal = matrix[j, j] - sum
+                    if (diagonal <= 0) {
+                        throw IllegalArgumentException("Matrix is not positive definite")
+                    }
+                    L[j, j] = sqrt(diagonal)
+                }
+                else -> {
+                    // Lower triangular elements
+                    var sum = 0.0
+                    for (k in 0 until j) {
+                        sum += L[i, k] * L[j, k]
+                    }
+                    L[i, j] = (matrix[i, j] - sum) / L[j, j]
+                }
+            }
+        }
+    }
+
+    return L
+}
+
+/**
+ * Solves the system of equations Ly = b using forward substitution, where L is a lower triangular matrix.
+ *
+ * @param L A lower triangular matrix of dimensions n x n.
+ * @param b A vector of size n representing the constant terms.
+ * @return A vector of size n representing the solution y to the system Ly = b.
+ */
+private fun forwardSubstitution(L: Matrix, b: DoubleArray): DoubleArray {
+    val n = L.rows
+    val y = DoubleArray(n)
+
+    for (i in 0 until n) {
+        var sum = 0.0
+        for (j in 0 until i) {
+            sum += L[i, j] * y[j]
+        }
+        y[i] = (b[i] - sum) / L[i, i]
+    }
+
+    return y
+}
+
+/**
+ * Performs backward substitution to solve a system of linear equations represented as L^T * x = y,
+ * where L is a lower triangular matrix and y is a known vector.
+ *
+ * This method assumes that the input matrix `L` is square and its diagonal contains no zeros.
+ *
+ * @param L A lower triangular matrix of type `Matrix` representing the system's coefficients.
+ *          It is accessed in a transposed fashion during the algorithm, where L[j, i] represents L^T[i, j].
+ * @param y A `DoubleArray` of size n representing the right-hand side vector of the linear system.
+ * @return A `DoubleArray` of size n representing the solution vector x such that L^T * x = y.
+ */
+private fun backwardSubstitution(L: Matrix, y: DoubleArray): DoubleArray {
+    val n = L.rows
+    val x = DoubleArray(n)
+
+    for (i in n - 1 downTo 0) {
+        var sum = 0.0
+        for (j in i + 1 until n) {
+            sum += L[j, i] * x[j]  // L[j, i] represents L^T[i, j]
+        }
+        x[i] = (y[i] - sum) / L[i, i]
+    }
+
+    return x
+}
+
+/**
+ * Computes the inverse of a given symmetric, positive-definite matrix using the Cholesky decomposition.
+ *
+ * This method assumes that the input matrix is square, symmetric, and positive-definite.
+ *
+ * @param matrix The input matrix to invert. Must be a square, symmetric, and positive-definite matrix.
+ * @return The inverse of the input matrix as a new `Matrix` object.
+ * @throws IllegalArgumentException if the input matrix is not square, symmetric, or positive-definite.
+ */
+fun invertMatrixCholesky(matrix: Matrix): Matrix {
+    val n = matrix.rows
+    val L = choleskyDecomposition(matrix)
+    val inverse = Matrix.zeros(n, n)
+
+    // For each column of the identity matrix
+    for (col in 0 until n) {
+        // Create unit vector for this column
+        val unitVector = DoubleArray(n)
+        unitVector[col] = 1.0
+
+        // Solve L * y = e_col
+        val y = forwardSubstitution(L, unitVector)
+
+        // Solve L^T * x = y
+        val x = backwardSubstitution(L, y)
+
+        // Store result in inverse matrix
+        for (row in 0 until n) {
+            inverse[row, col] = x[row]
+        }
+    }
+
+    return inverse
+}

orx-math/src/commonMain/kotlin/matrix/SparseMatrix.kt

@@ -0,0 +1,755 @@
+package org.openrndr.extra.math.matrix
+
+import kotlin.math.abs
+import kotlin.math.sqrt
+
+/**
+ * Counts the number of non-zero elements in a matrix.
+ *
+ * @param matrix The matrix to count non-zero elements in.
+ * @return The number of non-zero elements.
+ */
+private fun countNonZeroElements(matrix: Matrix): Int {
+    var count = 0
+    for (i in 0 until matrix.rows) {
+        for (j in 0 until matrix.cols) {
+            if (matrix[i, j] != 0.0) {
+                count++
+            }
+        }
+    }
+    return count
+}
+
+/**
+ * Represents a sparse matrix using the Compressed Sparse Row (CSR) format.
+ *
+ * The CSR format stores a sparse matrix using three arrays:
+ * - values: stores the non-zero values of the matrix
+ * - columnIndices: stores the column indices of the non-zero values
+ * - rowPointers: stores the starting position of each row in the values array
+ *
+ * This format is memory-efficient for sparse matrices where most elements are zero.
+ *
+ * @property rows The number of rows in the matrix.
+ * @property cols The number of columns in the matrix.
+ * @property values Array containing the non-zero values of the matrix.
+ * @property columnIndices Array containing the column indices of the non-zero values.
+ * @property rowPointers Array containing the starting position of each row in the values array.
+ */
+class SparseMatrix(
+    val rows: Int,
+    val cols: Int,
+    val values: DoubleArray,
+    val columnIndices: IntArray,
+    val rowPointers: IntArray
+) {
+
+
+    fun row(i: Int): DoubleArray {
+        val ret = DoubleArray(cols)
+        val start = rowPointers[i]
+        val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+        for (k in start until end) {
+            ret[columnIndices[k]] = values[k]
+        }
+
+        return ret
+    }
+
+
+    /**
+     * Gets the value at the specified position in the matrix.
+     *
+     * @param i The row index.
+     * @param j The column index.
+     * @return The value at position (i, j).
+     */
+    operator fun get(i: Int, j: Int): Double {
+        if (i < 0 || i >= rows || j < 0 || j >= cols) {
+            throw IndexOutOfBoundsException("Index out of bounds: ($i, $j)")
+        }
+
+        // Search for the value in the CSR format
+        val start = rowPointers[i]
+        val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+        for (k in start until end) {
+            if (columnIndices[k] == j) {
+                return values[k]
+            }
+        }
+
+        // If no value is found, return 0 (sparse matrix default)
+        return 0.0
+    }
+
+
+    /**
+     * Creates a dense Matrix representation of this sparse matrix.
+     *
+     * @return A dense Matrix with the same values as this sparse matrix.
+     */
+    fun toDenseMatrix(): Matrix {
+        val result = Matrix.zeros(rows, cols)
+
+        for (i in 0 until rows) {
+            val start = rowPointers[i]
+            val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+            for (k in start until end) {
+                val j = columnIndices[k]
+                result[i, j] = values[k]
+            }
+        }
+
+        return result
+    }
+
+
+    /**
+     * Multiplies this sparse matrix with another sparse matrix and returns the resulting sparse matrix.
+     *
+     * @param other The sparse matrix to multiply with this sparse matrix.
+     * @return The resulting sparse matrix after multiplication.
+     * @throws IllegalArgumentException If the number of columns in this matrix does not match the number of rows in the other matrix.
+     */
+    operator fun times(other: SparseMatrix): SparseMatrix {
+        if (cols != other.rows) {
+            throw IllegalArgumentException("Matrix dimensions don't match for multiplication")
+        }
+
+        // Create a temporary dense matrix to store the result
+        // We'll convert it to sparse at the end
+        val tempResult = Matrix.zeros(rows, other.cols)
+
+        // Perform the multiplication
+        for (i in 0 until rows) {
+            val rowStart = rowPointers[i]
+            val rowEnd = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+            for (k in rowStart until rowEnd) {
+                val col = columnIndices[k]
+                val value = values[k]
+
+                // For each non-zero element in the other matrix's row 'col'
+                val otherRowStart = other.rowPointers[col]
+                val otherRowEnd = if (col < other.rows - 1) other.rowPointers[col + 1] else other.values.size
+
+                for (l in otherRowStart until otherRowEnd) {
+                    val j = other.columnIndices[l]
+                    val otherValue = other.values[l]
+
+                    // Accumulate the product
+                    tempResult[i, j] += value * otherValue
+                }
+            }
+        }
+
+        // Convert the result to a sparse matrix
+        return fromMatrix(tempResult)
+    }
+
+
+    /**
+     * Multiplies this sparse matrix by a scalar value.
+     *
+     * @param scalar The scalar value to multiply by.
+     * @return A new sparse matrix with all elements multiplied by the scalar.
+     */
+    operator fun times(scalar: Double): SparseMatrix {
+        val newValues = DoubleArray(values.size) { i -> values[i] * scalar }
+        return SparseMatrix(rows, cols, newValues, columnIndices.copyOf(), rowPointers.copyOf())
+    }
+
+
+    /**
+     * Extension function to add two SparseMatrices.
+     *
+     * @param other The SparseMatrix to add.
+     * @return A new Matrix that is the sum of the two matrices.
+     * @throws IllegalArgumentException If the dimensions of the matrices don't match.
+     */
+    operator fun plus(other: SparseMatrix): SparseMatrix {
+        if (rows != other.rows || cols != other.cols) {
+            throw IllegalArgumentException("Matrix dimensions don't match for addition")
+        }
+
+        // Count non-zero elements in result
+        var nonZeroCount = 0
+        val countMap = mutableMapOf<Int, MutableMap<Int, Double>>()
+
+        // Process this matrix's non-zeros
+        for (i in 0 until rows) {
+            val start = rowPointers[i]
+            val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+            for (k in start until end) {
+                val j = columnIndices[k]
+                val value = values[k]
+                if (!countMap.containsKey(i)) {
+                    countMap[i] = mutableMapOf()
+                }
+                countMap[i]?.put(j, value)
+                nonZeroCount++
+            }
+        }
+
+        // Add other matrix's non-zeros 
+        for (i in 0 until other.rows) {
+            val start = other.rowPointers[i]
+            val end = if (i < other.rows - 1) other.rowPointers[i + 1] else other.values.size
+
+            for (k in start until end) {
+                val j = other.columnIndices[k]
+                val value = other.values[k]
+
+                if (countMap.containsKey(i) && countMap[i]?.containsKey(j) == true) {
+                    // Update existing value
+                    val newValue = countMap[i]?.get(j)!! + value
+                    countMap[i]?.put(j, newValue)
+                } else {
+                    // Add new value
+                    if (!countMap.containsKey(i)) {
+                        countMap[i] = mutableMapOf()
+                    }
+                    countMap[i]?.put(j, value)
+                    nonZeroCount++
+                }
+            }
+        }
+
+        // Create result arrays
+        val resultValues = DoubleArray(nonZeroCount)
+        val resultColIndices = IntArray(nonZeroCount)
+        val resultRowPointers = IntArray(rows + 1)
+
+        var valueIndex = 0
+        resultRowPointers[0] = 0
+
+        for (i in 0 until rows) {
+            val rowMap = countMap[i]
+            if (rowMap != null) {
+                for ((j, value) in rowMap.entries.sortedBy { it.key }) {
+                    resultValues[valueIndex] = value
+                    resultColIndices[valueIndex] = j
+                    valueIndex++
+                }
+            }
+            resultRowPointers[i + 1] = valueIndex
+        }
+
+        return SparseMatrix(rows, cols, resultValues, resultColIndices, resultRowPointers)
+    }
+
+    /**
+     * Extension function to subtract a SparseMatrix from another SparseMatrix.
+     *
+     * @param other The SparseMatrix to subtract.
+     * @return A new Matrix that is the difference of the two matrices.
+     * @throws IllegalArgumentException If the dimensions of the matrices don't match.
+     */
+    operator fun minus(other: SparseMatrix): SparseMatrix {
+        if (rows != other.rows || cols != other.cols) {
+            throw IllegalArgumentException("Matrix dimensions don't match for subtraction")
+        }
+
+
+        // Count non-zero elements in result
+        var nonZeroCount = 0
+        val countMap = mutableMapOf<Int, MutableMap<Int, Double>>()
+
+        // Process this matrix's non-zeros
+        for (i in 0 until rows) {
+            val start = rowPointers[i]
+            val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+            for (k in start until end) {
+                val j = columnIndices[k]
+                val value = values[k]
+                if (!countMap.containsKey(i)) {
+                    countMap[i] = mutableMapOf()
+                }
+                countMap[i]?.put(j, value)
+                nonZeroCount++
+            }
+        }
+
+        // Subtract other matrix's non-zeros
+        for (i in 0 until other.rows) {
+            val start = other.rowPointers[i]
+            val end = if (i < other.rows - 1) other.rowPointers[i + 1] else other.values.size
+
+            for (k in start until end) {
+                val j = other.columnIndices[k]
+                val value = other.values[k]
+
+                if (countMap.containsKey(i) && countMap[i]?.containsKey(j) == true) {
+                    // Update existing value
+                    val newValue = countMap[i]?.get(j)!! - value
+                    if (newValue != 0.0) {
+                        countMap[i]?.put(j, newValue)
+                    } else {
+                        countMap[i]?.remove(j)
+                        if (countMap[i]?.isEmpty() == true) {
+                            countMap.remove(i)
+                        }
+                        nonZeroCount--
+                    }
+                } else {
+                    // Add new negative value
+                    if (!countMap.containsKey(i)) {
+                        countMap[i] = mutableMapOf()
+                    }
+                    countMap[i]?.put(j, -value)
+                    nonZeroCount++
+                }
+            }
+        }
+
+        // Create result arrays
+        val resultValues = DoubleArray(nonZeroCount)
+        val resultColIndices = IntArray(nonZeroCount)
+        val resultRowPointers = IntArray(rows + 1)
+
+        var valueIndex = 0
+        resultRowPointers[0] = 0
+
+        for (i in 0 until rows) {
+            val rowMap = countMap[i]
+            if (rowMap != null) {
+                for ((j, value) in rowMap.entries.sortedBy { it.key }) {
+                    resultValues[valueIndex] = value
+                    resultColIndices[valueIndex] = j
+                    valueIndex++
+                }
+            }
+            resultRowPointers[i + 1] = valueIndex
+        }
+
+        return SparseMatrix(rows, cols, resultValues, resultColIndices, resultRowPointers)
+    }
+
+    /**
+     * Creates a copy of this sparse matrix.
+     *
+     * @return A new sparse matrix with the same values as this matrix.
+     */
+    fun copy(): SparseMatrix {
+        return SparseMatrix(
+            rows,
+            cols,
+            values.copyOf(),
+            columnIndices.copyOf(),
+            rowPointers.copyOf()
+        )
+    }
+
+    /**
+     * Checks if this sparse matrix is symmetric within a specified tolerance.
+     *
+     * A matrix is considered symmetric if it is square (same number of rows and columns)
+     * and satisfies the condition matrix[i, j] == matrix[j, i] for all i and j, within
+     * the given tolerance.
+     *
+     * @param tolerance The maximum allowable difference between matrix[i, j] and matrix[j, i]
+     *                  to still consider the matrix symmetric. Default is 1e-10.
+     * @return `true` if the matrix is symmetric within the specified tolerance, `false` otherwise.
+     */
+    fun isSymmetric(tolerance: Double = 1e-10): Boolean {
+        if (rows != cols) return false
+
+        // Convert to dense matrix for simplicity
+        // A more optimized implementation would check directly in the sparse format
+        val dense = toDenseMatrix()
+
+        for (i in 0 until rows) {
+            for (j in 0 until i) { // Only check lower triangle
+                if (abs(dense[i, j] - dense[j, i]) > tolerance) {
+                    return false
+                }
+            }
+        }
+        return true
+    }
+
+    /**
+     * Returns the number of non-zero elements in this sparse matrix.
+     *
+     * @return The number of non-zero elements.
+     */
+    fun nonZeroCount(): Int {
+        return values.size
+    }
+
+    /**
+     * Returns the density of this sparse matrix (ratio of non-zero elements to total elements).
+     *
+     * @return The density as a value between 0 and 1.
+     */
+    fun density(): Double {
+        return values.size.toDouble() / (rows * cols)
+    }
+
+    /**
+     * Creates a transposed version of this sparse matrix.
+     *
+     * @return A new SparseMatrix that is the transpose of this matrix.
+     */
+    fun transpose(): SparseMatrix {
+        // Count non-zero elements per column for creating row pointers in transposed matrix
+        val colCounts = IntArray(cols)
+        for (k in columnIndices.indices) {
+            colCounts[columnIndices[k]]++
+        }
+
+        // Create row pointers for transposed matrix
+        val newRowPointers = IntArray(cols + 1)
+        newRowPointers[0] = 0
+        for (j in 0 until cols) {
+            newRowPointers[j + 1] = newRowPointers[j] + colCounts[j]
+        }
+
+        val newValues = DoubleArray(values.size)
+        val newColumnIndices = IntArray(values.size)
+
+        // Keep track of current position in each column
+        val colPositions = newRowPointers.copyOf()
+
+        // Fill in values and column indices for transposed matrix
+        for (i in 0 until rows) {
+            val start = rowPointers[i]
+            val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+            for (k in start until end) {
+                val j = columnIndices[k]
+                val pos = colPositions[j]
+                newValues[pos] = values[k]
+                newColumnIndices[pos] = i
+                colPositions[j]++
+            }
+        }
+
+        return SparseMatrix(cols, rows, newValues, newColumnIndices, newRowPointers)
+    }
+
+    companion object {
+        /**
+         * Creates a sparse matrix from a dense matrix.
+         *
+         * @param matrix The dense matrix to convert.
+         * @return A sparse matrix representation of the input matrix.
+         */
+        fun fromMatrix(matrix: Matrix): SparseMatrix {
+            val nonZeroCount = countNonZeroElements(matrix)
+
+            val values = DoubleArray(nonZeroCount)
+            val columnIndices = IntArray(nonZeroCount)
+            val rowPointers = IntArray(matrix.rows + 1)
+
+            var valueIndex = 0
+            rowPointers[0] = 0
+
+            for (i in 0 until matrix.rows) {
+                for (j in 0 until matrix.cols) {
+                    val value = matrix[i, j]
+                    if (value != 0.0) {
+                        values[valueIndex] = value
+                        columnIndices[valueIndex] = j
+                        valueIndex++
+                    }
+                }
+                rowPointers[i + 1] = valueIndex
+            }
+
+            return SparseMatrix(matrix.rows, matrix.cols, values, columnIndices, rowPointers)
+        }
+
+        /**
+         * Creates a sparse identity matrix of the specified size.
+         *
+         * @param n The size of the identity matrix (number of rows and columns).
+         * @return A sparse identity matrix of size n x n.
+         */
+        fun identity(n: Int): SparseMatrix {
+            val values = DoubleArray(n) { 1.0 }
+            val columnIndices = IntArray(n) { it }
+            val rowPointers = IntArray(n + 1) { it }
+
+            return SparseMatrix(n, n, values, columnIndices, rowPointers)
+        }
+
+        /**
+         * Creates a sparse zero matrix with the specified dimensions.
+         *
+         * @param rows The number of rows in the matrix.
+         * @param cols The number of columns in the matrix.
+         * @return A sparse zero matrix with dimensions [rows] x [cols].
+         */
+        fun zeros(rows: Int, cols: Int): SparseMatrix {
+            // For a zero matrix, there are no non-zero elements
+            val values = DoubleArray(0)
+            val columnIndices = IntArray(0)
+            val rowPointers = IntArray(rows + 1) { 0 }
+
+            return SparseMatrix(rows, cols, values, columnIndices, rowPointers)
+        }
+    }
+}
+
+/**
+ * Extension function to convert a dense Matrix to a SparseMatrix.
+ *
+ * @return A sparse matrix representation of this dense matrix.
+ */
+fun Matrix.toSparseMatrix(): SparseMatrix {
+    return SparseMatrix.fromMatrix(this)
+}
+
+
+/**
+ * Extension function to multiply a scalar by a SparseMatrix.
+ *
+ * @param matrix The SparseMatrix to multiply.
+ * @return A new SparseMatrix with all elements multiplied by the scalar.
+ */
+operator fun Double.times(matrix: SparseMatrix): SparseMatrix {
+    return matrix * this
+}
+
+
+/**
+ * Extension function to add a Matrix to a SparseMatrix.
+ *
+ * @param other The SparseMatrix to add.
+ * @return A new Matrix that is the sum of the two matrices.
+ * @throws IllegalArgumentException If the dimensions of the matrices don't match.
+ */
+operator fun Matrix.plus(other: SparseMatrix): Matrix {
+    if (rows != other.rows || cols != other.cols) {
+        throw IllegalArgumentException("Matrix dimensions don't match for addition")
+    }
+
+    val result = this.copy()
+    val otherDense = other.toDenseMatrix()
+
+    for (i in 0 until rows) {
+        for (j in 0 until cols) {
+            result[i, j] += otherDense[i, j]
+        }
+    }
+
+    return result
+}
+
+/**
+ * Extension function to subtract a Matrix from a SparseMatrix.
+ *
+ * @param other The Matrix to subtract.
+ * @return A new Matrix that is the difference of the two matrices.
+ * @throws IllegalArgumentException If the dimensions of the matrices don't match.
+ */
+operator fun SparseMatrix.minus(other: Matrix): Matrix {
+    if (rows != other.rows || cols != other.cols) {
+        throw IllegalArgumentException("Matrix dimensions don't match for subtraction")
+    }
+
+    val result = toDenseMatrix()
+
+    for (i in 0 until rows) {
+        for (j in 0 until cols) {
+            result[i, j] -= other[i, j]
+        }
+    }
+
+    return result
+}
+
+/**
+ * Adds this sparse matrix to another matrix and returns the resulting matrix.
+ *
+ * @param other The matrix to add to this sparse matrix.
+ * @return The resulting matrix after addition.
+ * @throws IllegalArgumentException If the dimensions of the matrices don't match.
+ */
+operator fun SparseMatrix.plus(other: Matrix): Matrix {
+    if (rows != other.rows || cols != other.cols) {
+        throw IllegalArgumentException("Matrix dimensions don't match for addition")
+    }
+
+    // Convert to dense matrix for simplicity
+    val result = toDenseMatrix()
+
+    for (i in 0 until rows) {
+        for (j in 0 until cols) {
+            result[i, j] += other[i, j]
+        }
+    }
+
+    return result
+}
+
+/**
+ * Extension function to subtract a SparseMatrix from a Matrix.
+ *
+ * @param other The SparseMatrix to subtract.
+ * @return A new Matrix that is the difference of the two matrices.
+ * @throws IllegalArgumentException If the dimensions of the matrices don't match.
+ */
+operator fun Matrix.minus(other: SparseMatrix): Matrix {
+    if (rows != other.rows || cols != other.cols) {
+        throw IllegalArgumentException("Matrix dimensions don't match for subtraction")
+    }
+
+    val result = this.copy()
+    val otherDense = other.toDenseMatrix()
+
+    for (i in 0 until rows) {
+        for (j in 0 until cols) {
+            result[i, j] -= otherDense[i, j]
+        }
+    }
+
+    return result
+}
+
+
+/**
+ * Multiplies this sparse matrix with another matrix and returns the resulting matrix.
+ *
+ * @param other The matrix to multiply with this sparse matrix.
+ * @return The resulting matrix after multiplication.
+ * @throws IllegalArgumentException If the number of columns in this matrix does not match the number of rows in the other matrix.
+ */
+operator fun SparseMatrix.times(other: Matrix): Matrix {
+    if (cols != other.rows) {
+        throw IllegalArgumentException("Matrix dimensions don't match for multiplication")
+    }
+
+    val result = Matrix.zeros(rows, other.cols)
+
+    for (i in 0 until rows) {
+        val start = rowPointers[i]
+        val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+        for (j in 0 until other.cols) {
+            var sum = 0.0
+            for (k in start until end) {
+                val col = columnIndices[k]
+                sum += values[k] * other[col, j]
+            }
+            result[i, j] = sum
+        }
+    }
+
+    return result
+}
+
+/**
+ * Constructs a SparseMatrix instance with the specified number of rows, columns, and non-zero values.
+ *
+ * @param rows The number of rows in the matrix.
+ * @param cols The number of columns in the matrix.
+ * @param values A map containing the non-zero values and their positions, where the key is a pair of row and column indices, and the value is the corresponding matrix value.
+ * @return A new SparseMatrix instance with the given dimensions and non-zero values.
+ */
+fun SparseMatrix(rows: Int, cols: Int, values: Map<Pair<Int, Int>, Double>): SparseMatrix {
+    val indices = values.entries.map { it.key }
+    val values = values.entries.map { it.value }
+    return SparseMatrix(rows, cols, indices, values, true)
+}
+
+/**
+ * Constructs a sparse matrix in Compressed Sparse Row (CSR) format.
+ *
+ * @param rows The number of rows in the matrix.
+ * @param cols The number of columns in the matrix.
+ * @param indices A list of pairs representing the row and column indices of non-zero elements.
+ * @param values A list of values corresponding to the non-zero elements at the specified indices.
+ * @param presorted Indicates whether the indices are already sorted by row and column. If set to `false`, the indices
+ *        will be sorted internally. Defaults to `false`.
+ * @return The constructed `SparseMatrix` object in CSR format.
+ * @throws IllegalArgumentException If the number of indices does not match the number of values, or if the indices list is empty.
+ */
+
+fun SparseMatrix(
+    rows: Int,
+    cols: Int,
+    indices: List<Pair<Int, Int>>,
+    values: List<Double>,
+    presorted: Boolean = true
+): SparseMatrix {
+    if (indices.size != values.size) {
+        throw IllegalArgumentException("Number of indices must match number of values")
+    }
+    if (indices.isEmpty()) {
+        throw IllegalArgumentException("Cannot create sparse matrix from empty data")
+    }
+
+    // Sort entries by row, then column
+    val sortedEntries = indices.zip(values).let {
+        if (presorted) it else it.sortedWith(compareBy({ it.first.first }, { it.first.second }))
+    }
+
+    // Create arrays for CSR format
+    val resultValues = DoubleArray(values.size)
+    val resultColIndices = IntArray(values.size)
+    val resultRowPointers = IntArray(rows + 1)
+
+    // Fill in values and column indices
+    for (i in sortedEntries.indices) {
+        val entry = sortedEntries[i]
+        resultValues[i] = entry.second
+        resultColIndices[i] = entry.first.second
+    }
+
+    // Create row pointers
+    var currentRow = 0
+    resultRowPointers[0] = 0
+    for (i in sortedEntries.indices) {
+        val row = sortedEntries[i].first.first
+        while (currentRow < row) {
+            currentRow++
+            resultRowPointers[currentRow] = i
+        }
+    }
+    while (currentRow < rows) {
+        currentRow++
+        resultRowPointers[currentRow] = sortedEntries.size
+    }
+
+    return SparseMatrix(rows, cols, resultValues, resultColIndices, resultRowPointers)
+}
+
+fun SparseMatrix.checkIntegrity() : SparseMatrix {
+    for (i in 0 until rows) {
+        val start = rowPointers[i]
+        val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+
+        val indices = (start until end).map {
+            columnIndices[it]
+        }
+        require(indices.size == indices.distinct().size) { "Duplicate indices found in row $i" }
+    }
+    return this
+}
+
+/**
+ * Calculates the Frobenius norm of the sparse matrix.
+ * The Frobenius norm is the square root of the sum of the squares
+ * of all the non-zero elements in the matrix.
+ *
+ * @return The Frobenius norm as a Double value.
+ * @see
+ */
+fun SparseMatrix.frobeniusNorm() : Double {
+    var sum = 0.0
+    for (i in 0 until rows) {
+        val start = rowPointers[i]
+        val end = if (i < rows - 1) rowPointers[i + 1] else values.size
+        for (j in start until end) {
+            sum += values[j] * values[j]
+        }
+    }
+    return sqrt(sum)
+}

orx-math/src/commonMain/kotlin/matrix/SparseMatrixQR.kt

@@ -0,0 +1,148 @@
+package org.openrndr.extra.math.matrix
+
+import kotlin.math.abs
+import kotlin.math.sign
+import kotlin.math.sqrt
+
+
+/**
+ * Performs QR decomposition on a sparse matrix.
+ *
+ * QR decomposition factors a matrix A into the product of an orthogonal matrix Q
+ * and an upper triangular matrix R, such that A = QR.
+ *
+ * This implementation uses the Gram-Schmidt process adapted for sparse matrices.
+ * It's more straightforward than Householder reflections and easier to debug.
+ *
+ * This implementation is optimized for sparse matrices by:
+ * 1. Efficiently extracting columns from the sparse matrix
+ * 2. Avoiding unnecessary operations on zero elements
+ * 3. Returning sparse matrices as the result
+ *
+ * @param matrix The sparse matrix to decompose.
+ * @return A Pair of matrices (Q, R) representing the orthogonal and upper triangular matrices.
+ */
+fun qrDecomposition(matrix: SparseMatrix): Pair<SparseMatrix, SparseMatrix> {
+    val m = matrix.rows
+    val n = matrix.cols
+
+    // Initialize Q and R matrices
+    val q = Matrix(m, n)
+    val r = Matrix(n, n)
+
+    // Implement the Gram-Schmidt process
+    for (j in 0 until n) {
+        // Extract the j-th column of the sparse matrix efficiently
+        for (i in 0 until m) {
+            q[i, j] = matrix[i, j]
+        }
+
+        // Orthogonalize against previous columns
+        for (k in 0 until j) {
+            // Compute dot product of q_j and q_k
+            var dotProduct = 0.0
+            for (i in 0 until m) {
+                // Only multiply non-zero elements
+                if (q[i, j] != 0.0 && q[i, k] != 0.0) {
+                    dotProduct += q[i, j] * q[i, k]
+                }
+            }
+
+            // Store in R
+            r[k, j] = dotProduct
+
+            // Subtract projection (only for non-zero elements in q_k)
+            for (i in 0 until m) {
+                if (q[i, k] != 0.0) {
+                    q[i, j] -= dotProduct * q[i, k]
+                }
+            }
+        }
+
+        // Compute the norm of the column (only for non-zero elements)
+        var norm = 0.0
+        for (i in 0 until m) {
+            if (q[i, j] != 0.0) {
+                norm += q[i, j] * q[i, j]
+            }
+        }
+        norm = sqrt(norm)
+
+        // Store the norm in R
+        r[j, j] = norm
+
+        // Normalize the column if it's not too small
+        if (norm > 1e-10) {
+            for (i in 0 until m) {
+                if (q[i, j] != 0.0) {
+                    q[i, j] /= norm
+                }
+            }
+        } else {
+            // If the column is too small, set it to zero
+            for (i in 0 until m) {
+                q[i, j] = 0.0
+            }
+        }
+    }
+
+    // Convert matrices back to sparse format
+    val qSparse = q.toSparseMatrix()
+    val rSparse = r.toSparseMatrix()
+
+    return Pair(qSparse, rSparse)
+}
+
+/**
+ * Solves a linear system Ax = b using QR decomposition.
+ *
+ * Given the QR decomposition of matrix A (where A = QR), this function solves
+ * the system Ax = b by first computing y = Q^T * b, and then solving Rx = y
+ * using backward substitution.
+ *
+ * @param qr A Pair of sparse matrices (Q, R) representing the QR decomposition of matrix A.
+ * @param b The right-hand side vector/matrix of the system.
+ * @return The solution vector/matrix x.
+ * @throws IllegalArgumentException If the dimensions of the matrices are incompatible.
+ */
+fun solveQR(qr: Pair<SparseMatrix, SparseMatrix>, b: SparseMatrix): SparseMatrix {
+    val (q, r) = qr
+
+    // Check dimensions
+    if (q.rows != b.rows) {
+        throw IllegalArgumentException("Incompatible dimensions: Q has ${q.rows} rows, b has ${b.rows} rows")
+    }
+
+    val n = r.rows  // Number of columns in the original matrix
+    val m = b.cols  // Number of right-hand sides
+
+    // Compute y = Q^T * b
+    val qTranspose = q.transpose()
+    val y = qTranspose * b
+
+    // Create a dense matrix to store the solution
+    val x = Matrix.zeros(n, m)
+
+    // Backward substitution: Solve Rx = y for x
+    for (col in 0 until m) {  // For each right-hand side
+        for (i in n - 1 downTo 0) {  // Start from the last row
+            var sum = y[i, col]
+
+            // Subtract the known terms
+            for (k in i + 1 until n) {
+                sum -= r[i, k] * x[k, col]
+            }
+
+            // Divide by the diagonal element of R
+            val rii = r[i, i]
+            if (abs(rii) < 1e-10) {
+                throw IllegalArgumentException("R matrix is singular or nearly singular")
+            }
+
+            x[i, col] = sum / rii
+        }
+    }
+
+    // Convert to sparse matrix
+    return x.toSparseMatrix()
+}

orx-math/src/commonTest/kotlin/matrix/SparseMatrixQRTest.kt

@@ -0,0 +1,123 @@
+package matrix
+
+import org.openrndr.extra.math.matrix.SparseMatrix
+import org.openrndr.extra.math.matrix.qrDecomposition
+import org.openrndr.extra.math.matrix.solveQR
+import kotlin.test.assertEquals
+import kotlin.math.abs
+import kotlin.test.Test
+
+class SparseMatrixQRTest {
+
+    @Test
+    fun testQRDecomposition() {
+        // Create a sparse matrix
+        // This will create the following matrix:
+        // 4.0, 3.0, 0.0
+        // 6.0, 3.0, 0.0
+        // 0.0, 1.0, 5.0
+        val values = doubleArrayOf(4.0, 3.0, 6.0, 3.0, 1.0, 5.0)
+        val columnIndices = intArrayOf(0, 1, 0, 1, 1, 2)
+        val rowPointers = intArrayOf(0, 2, 4, 6)
+        val sparseMatrix = SparseMatrix(3, 3, values, columnIndices, rowPointers)
+
+        // Perform QR decomposition
+        val (q, r) = qrDecomposition(sparseMatrix)
+
+        val originalDense = sparseMatrix.toDenseMatrix()
+        for (i in 0 until originalDense.rows) {
+            val row = (0 until originalDense.cols).map { j -> originalDense[i, j] }.joinToString()
+        }
+
+        val qDense = q.toDenseMatrix()
+        for (i in 0 until qDense.rows) {
+            val row = (0 until qDense.cols).map { j -> qDense[i, j] }.joinToString()
+        }
+
+        val rDense = r.toDenseMatrix()
+        for (i in 0 until rDense.rows) {
+            val row = (0 until rDense.cols).map { j -> rDense[i, j] }.joinToString()
+        }
+
+        // Verify Q is orthogonal (Q^T * Q = I)
+        val qTranspose = q.transpose()
+        val qTq = qTranspose * q
+        val identity = SparseMatrix.identity(q.rows)
+
+        val qTqDense = qTq.toDenseMatrix()
+        for (i in 0 until qTqDense.rows) {
+            val row = (0 until qTqDense.cols).map { j -> qTqDense[i, j] }.joinToString()
+        }
+
+        // Check that Q^T * Q is approximately identity
+        for (i in 0 until q.rows) {
+            for (j in 0 until q.cols) {
+                assertEquals(identity[i, j], qTq[i, j], 1e-10)
+            }
+        }
+
+        // Verify R is upper triangular
+        for (i in 0 until r.rows) {
+            for (j in 0 until r.cols) {
+                if (i > j) {
+                    assertEquals(0.0, abs(r[i, j]), 1e-10)
+                }
+            }
+        }
+
+        // Verify A = Q * R
+        val qr = q * r
+        
+        val qrDense = qr.toDenseMatrix()
+        for (i in 0 until qrDense.rows) {
+            val row = (0 until qrDense.cols).map { j -> qrDense[i, j] }.joinToString()
+        }
+
+        // Check that A = Q * R
+        for (i in 0 until sparseMatrix.rows) {
+            for (j in 0 until sparseMatrix.cols) {
+                assertEquals(sparseMatrix[i, j], qr[i, j], 1e-10)
+            }
+        }
+    }
+
+    @Test
+    fun testSolveWithQRDecomposition() {
+        // Create a sparse matrix A
+        // This will create the following matrix:
+        // 4.0, 3.0, 0.0
+        // 6.0, 3.0, 0.0
+        // 0.0, 1.0, 5.0
+        val valuesA = doubleArrayOf(4.0, 3.0, 6.0, 3.0, 1.0, 5.0)
+        val columnIndicesA = intArrayOf(0, 1, 0, 1, 1, 2)
+        val rowPointersA = intArrayOf(0, 2, 4, 6)
+        val matrixA = SparseMatrix(3, 3, valuesA, columnIndicesA, rowPointersA)
+
+        // Create a sparse matrix b (right-hand side)
+        // This will create the following vector:
+        // 1.0
+        // 2.0
+        // 3.0
+        val valuesB = doubleArrayOf(1.0, 2.0, 3.0)
+        val columnIndicesB = intArrayOf(0, 0, 0)
+        val rowPointersB = intArrayOf(0, 1, 2, 3)
+        val matrixB = SparseMatrix(3, 1, valuesB, columnIndicesB, rowPointersB)
+
+        // Perform QR decomposition
+        val qr = qrDecomposition(matrixA)
+
+        // Solve the system Ax = b
+        val x = solveQR(qr, matrixB)
+
+
+        // Verify A * x = b
+        val product = matrixA * x
+
+        // Check that A * x = b
+        for (i in 0 until matrixB.rows) {
+            for (j in 0 until matrixB.cols) {
+                assertEquals(matrixB[i, j], product[i, j], 1e-10)
+            }
+        }
+    }
+}

orx-math/src/commonTest/kotlin/matrix/SparseMatrixTest.kt

@@ -0,0 +1,270 @@
+package matrix
+
+import org.openrndr.extra.math.matrix.Matrix
+import org.openrndr.extra.math.matrix.SparseMatrix
+import org.openrndr.extra.math.matrix.checkIntegrity
+import org.openrndr.extra.math.matrix.toSparseMatrix
+import kotlin.test.Test
+import kotlin.test.assertEquals
+import kotlin.test.assertFalse
+import kotlin.test.assertTrue
+
+class SparseMatrixTest {
+
+    @Test
+    fun testCreateSparseMatrix() {
+        // Create a sparse matrix directly
+        val values = doubleArrayOf(1.0, 2.0, 3.0, 4.0)
+        val columnIndices = intArrayOf(0, 2, 1, 2)
+        val rowPointers = intArrayOf(0, 2, 3, 4)
+        val sparseMatrix = SparseMatrix(3, 3, values, columnIndices, rowPointers)
+
+        // Verify dimensions
+        assertEquals(3, sparseMatrix.rows)
+        assertEquals(3, sparseMatrix.cols)
+
+        // Verify values
+        assertEquals(1.0, sparseMatrix[0, 0])
+        assertEquals(0.0, sparseMatrix[0, 1])
+        assertEquals(2.0, sparseMatrix[0, 2])
+        assertEquals(0.0, sparseMatrix[1, 0])
+        assertEquals(3.0, sparseMatrix[1, 1])
+        assertEquals(0.0, sparseMatrix[1, 2])
+        assertEquals(0.0, sparseMatrix[2, 0])
+        assertEquals(0.0, sparseMatrix[2, 1])
+        assertEquals(4.0, sparseMatrix[2, 2])
+    }
+
+    @Test
+    fun testFromDenseMatrix() {
+        // Create a dense matrix
+        val denseMatrix = Matrix(3, 3)
+        denseMatrix[0, 0] = 1.0
+        denseMatrix[0, 2] = 2.0
+        denseMatrix[1, 1] = 3.0
+        denseMatrix[2, 2] = 4.0
+
+        // Convert to sparse matrix
+        val sparseMatrix = denseMatrix.toSparseMatrix()
+
+        // Verify dimensions
+        assertEquals(3, sparseMatrix.rows)
+        assertEquals(3, sparseMatrix.cols)
+
+        // Verify values
+        assertEquals(1.0, sparseMatrix[0, 0])
+        assertEquals(0.0, sparseMatrix[0, 1])
+        assertEquals(2.0, sparseMatrix[0, 2])
+        assertEquals(0.0, sparseMatrix[1, 0])
+        assertEquals(3.0, sparseMatrix[1, 1])
+        assertEquals(0.0, sparseMatrix[1, 2])
+        assertEquals(0.0, sparseMatrix[2, 0])
+        assertEquals(0.0, sparseMatrix[2, 1])
+        assertEquals(4.0, sparseMatrix[2, 2])
+
+        // Verify non-zero count
+        assertEquals(4, sparseMatrix.nonZeroCount())
+    }
+
+    @Test
+    fun testToDenseMatrix() {
+        // Create a sparse matrix
+        val values = doubleArrayOf(1.0, 2.0, 3.0, 4.0)
+        val columnIndices = intArrayOf(0, 2, 1, 2)
+        val rowPointers = intArrayOf(0, 2, 3, 4)
+        val sparseMatrix = SparseMatrix(3, 3, values, columnIndices, rowPointers)
+
+        // Convert to dense matrix
+        val denseMatrix = sparseMatrix.toDenseMatrix()
+
+        // Verify dimensions
+        assertEquals(3, denseMatrix.rows)
+        assertEquals(3, denseMatrix.cols)
+
+        // Verify values
+        assertEquals(1.0, denseMatrix[0, 0])
+        assertEquals(0.0, denseMatrix[0, 1])
+        assertEquals(2.0, denseMatrix[0, 2])
+        assertEquals(0.0, denseMatrix[1, 0])
+        assertEquals(3.0, denseMatrix[1, 1])
+        assertEquals(0.0, denseMatrix[1, 2])
+        assertEquals(0.0, denseMatrix[2, 0])
+        assertEquals(0.0, denseMatrix[2, 1])
+        assertEquals(4.0, denseMatrix[2, 2])
+    }
+
+    @Test
+    fun testMatrixMultiplication() {
+        // Create a sparse matrix
+        val values1 = doubleArrayOf(1.0, 2.0, 3.0)
+        val columnIndices1 = intArrayOf(0, 1, 0)
+        val rowPointers1 = intArrayOf(0, 2, 3)
+        val sparseMatrix1 = SparseMatrix(2, 2, values1, columnIndices1, rowPointers1)
+
+        // Create another sparse matrix
+        val values2 = doubleArrayOf(4.0, 5.0, 6.0)
+        val columnIndices2 = intArrayOf(0, 1, 1)
+        val rowPointers2 = intArrayOf(0, 2, 3)
+        val sparseMatrix2 = SparseMatrix(2, 2, values2, columnIndices2, rowPointers2)
+
+        // Multiply sparse matrices
+        val result = sparseMatrix1.times(sparseMatrix2)
+
+        result.checkIntegrity()
+
+        // Verify dimensions
+        assertEquals(2, result.rows)
+        assertEquals(2, result.cols)
+
+        // Verify values (1*4 + 2*0 = 4, 1*5 + 2*6 = 17, 3*4 + 0*0 = 12, 3*5 + 0*6 = 15)
+        assertEquals(4.0, result[0, 0])
+        assertEquals(17.0, result[0, 1])
+        assertEquals(12.0, result[1, 0])
+        assertEquals(15.0, result[1, 1])
+    }
+
+    @Test
+    fun testScalarMultiplication() {
+        // Create a sparse matrix
+        val values = doubleArrayOf(1.0, 2.0, 3.0, 4.0)
+        val columnIndices = intArrayOf(0, 2, 1, 2)
+        val rowPointers = intArrayOf(0, 2, 3, 4)
+        val sparseMatrix = SparseMatrix(3, 3, values, columnIndices, rowPointers)
+
+        // Multiply by scalar
+        val result = sparseMatrix * 2.0
+
+        // Verify dimensions
+        assertEquals(3, result.rows)
+        assertEquals(3, result.cols)
+
+        // Verify values
+        assertEquals(2.0, result[0, 0])
+        assertEquals(0.0, result[0, 1])
+        assertEquals(4.0, result[0, 2])
+        assertEquals(0.0, result[1, 0])
+        assertEquals(6.0, result[1, 1])
+        assertEquals(0.0, result[1, 2])
+        assertEquals(0.0, result[2, 0])
+        assertEquals(0.0, result[2, 1])
+        assertEquals(8.0, result[2, 2])
+    }
+
+    @Test
+    fun testAddition() {
+        // Create a sparse matrix
+        val values1 = doubleArrayOf(1.0, 2.0, 3.0)
+        val columnIndices1 = intArrayOf(0, 2, 1)
+        val rowPointers1 = intArrayOf(0, 2, 3)
+        val sparseMatrix1 = SparseMatrix(2, 3, values1, columnIndices1, rowPointers1)
+
+        // Create another sparse matrix
+        val values2 = doubleArrayOf(4.0, 5.0, 6.0)
+        val columnIndices2 = intArrayOf(0, 1, 2)
+        val rowPointers2 = intArrayOf(0, 2, 3)
+        val sparseMatrix2 = SparseMatrix(2, 3, values2, columnIndices2, rowPointers2)
+
+        // Add sparse matrices
+        val result = sparseMatrix1 + sparseMatrix2
+        result.checkIntegrity()
+
+        // Verify dimensions
+        assertEquals(2, result.rows)
+        assertEquals(3, result.cols)
+
+        // Verify values
+        assertEquals(5.0, result[0, 0])
+        assertEquals(5.0, result[0, 1])
+        assertEquals(2.0, result[0, 2])
+        assertEquals(0.0, result[1, 0])
+        assertEquals(3.0, result[1, 1])
+        assertEquals(6.0, result[1, 2])
+    }
+
+    @Test
+    fun testSubtraction() {
+        // Create a sparse matrix
+        val values1 = doubleArrayOf(5.0, 7.0, 9.0)
+        val columnIndices1 = intArrayOf(0, 2, 1)
+        val rowPointers1 = intArrayOf(0, 2, 3)
+        val sparseMatrix1 = SparseMatrix(2, 3, values1, columnIndices1, rowPointers1)
+
+        // Create another sparse matrix
+        val values2 = doubleArrayOf(1.0, 2.0, 3.0)
+        val columnIndices2 = intArrayOf(0, 1, 2)
+        val rowPointers2 = intArrayOf(0, 2, 3)
+        val sparseMatrix2 = SparseMatrix(2, 3, values2, columnIndices2, rowPointers2)
+
+        // Subtract sparse matrices
+        val result = sparseMatrix1 - sparseMatrix2
+        result.checkIntegrity()
+
+        // Verify dimensions
+        assertEquals(2, result.rows)
+        assertEquals(3, result.cols)
+
+        // Verify values
+        assertEquals(4.0, result[0, 0])
+        assertEquals(-2.0, result[0, 1])
+        assertEquals(7.0, result[0, 2])
+        assertEquals(0.0, result[1, 0])
+        assertEquals(9.0, result[1, 1])
+        assertEquals(-3.0, result[1, 2])
+    }
+
+    @Test
+    fun testTranspose() {
+        // Create a sparse matrix
+        val values = doubleArrayOf(1.0, 2.0, 3.0, 4.0)
+        val columnIndices = intArrayOf(0, 2, 1, 2)
+        val rowPointers = intArrayOf(0, 2, 3, 4)
+        val sparseMatrix = SparseMatrix(3, 3, values, columnIndices, rowPointers)
+        sparseMatrix.checkIntegrity()
+
+        // Transpose the matrix
+        val transposed = sparseMatrix.transpose()
+
+        // Verify dimensions
+        assertEquals(3, transposed.rows)
+        assertEquals(3, transposed.cols)
+
+        // Verify values
+        assertEquals(1.0, transposed[0, 0])
+        assertEquals(0.0, transposed[0, 1])
+        assertEquals(0.0, transposed[0, 2])
+        assertEquals(0.0, transposed[1, 0])
+        assertEquals(3.0, transposed[1, 1])
+        assertEquals(0.0, transposed[1, 2])
+        assertEquals(2.0, transposed[2, 0])
+        assertEquals(0.0, transposed[2, 1])
+        assertEquals(4.0, transposed[2, 2])
+    }
+
+    @Test
+    fun testIsSymmetric() {
+        // Create a symmetric sparse matrix
+        // For a symmetric matrix, if matrix[i,j] = value, then matrix[j,i] = value
+        // In CSR format, we need to ensure both positions have the same value
+
+        // This will create the following symmetric matrix:
+        // 1.0, 2.0, 3.0
+        // 2.0, 4.0, 5.0
+        // 3.0, 5.0, 6.0
+        val values = doubleArrayOf(1.0, 2.0, 3.0, 2.0, 4.0, 5.0, 3.0, 5.0, 6.0)
+        val columnIndices = intArrayOf(0, 1, 2, 0, 1, 2, 0, 1, 2)
+        val rowPointers = intArrayOf(0, 3, 6)
+        val symmetricMatrix = SparseMatrix(3, 3, values, columnIndices, rowPointers)
+
+        // Verify it's symmetric
+        assertTrue(symmetricMatrix.isSymmetric())
+
+        // Create a non-symmetric sparse matrix
+        val values2 = doubleArrayOf(1.0, 2.0, 3.0, 4.0)
+        val columnIndices2 = intArrayOf(0, 2, 1, 2)
+        val rowPointers2 = intArrayOf(0, 2, 3, 4)
+        val nonSymmetricMatrix = SparseMatrix(3, 3, values2, columnIndices2, rowPointers2)
+
+        // Verify it's not symmetric
+        assertFalse(nonSymmetricMatrix.isSymmetric())
+    }
+}

orx-math/src/jvmDemo/kotlin/matrix/DemoLeastSquares01.kt

@@ -0,0 +1,65 @@
+package matrix
+
+import org.openrndr.application
+import org.openrndr.color.ColorRGBa
+import org.openrndr.extra.math.matrix.Matrix
+import org.openrndr.extra.math.matrix.invertMatrixCholesky
+import org.openrndr.extra.noise.uniform
+import org.openrndr.extra.noise.uniformRing
+import org.openrndr.math.Vector2
+import kotlin.math.cos
+import kotlin.random.Random
+
+/**
+ * Demonstrate least squares method to find a regression line through noisy points
+ * Line drawn in red is the estimated line, in green is the ground-truth line
+ *
+ * Ax = b => x = A⁻¹b
+ * because A is likely inconsistent, we look for an approximate x based on AᵀA, which is consistent.
+ * x̂ = (AᵀA)⁻¹ Aᵀb
+ */
+fun main() {
+    application {
+        configure {
+            width = 720
+            height = 720
+        }
+        program {
+            extend {
+                val ls = drawer.bounds.horizontal(0.5).rotateBy(cos(seconds)*45.0)
+
+                val r = Random((seconds*10).toInt())
+
+                val pointCount = 100
+                val A = Matrix(pointCount, 2)
+                val b = Matrix(pointCount, 1)
+                for (i in 0 until pointCount) {
+
+                    val p = ls.position(Double.uniform(0.0, 1.0, r))
+                    val pr = p + Vector2.uniformRing(0.0, 130.0, r)
+
+                    A[i, 0] = 1.0
+                    A[i, 1] = pr.x
+                    b[i, 0] = pr.y
+
+                    drawer.circle(pr, 5.0)
+                }
+                val At = A.transposed()
+                val AtA = At * A
+                val Atb = At * b
+
+                val AtAI = invertMatrixCholesky(AtA)
+                val x = AtAI * Atb
+
+                val p0 = Vector2(0.0, x[0, 0])
+                val p1 = Vector2(720.0, x[0, 0] + x[1, 0] * 720.0)
+
+                drawer.stroke = ColorRGBa.RED
+                drawer.lineSegment(p0, p1)
+
+                drawer.stroke = ColorRGBa.GREEN
+                drawer.lineSegment(ls)
+            }
+        }
+    }
+}