KMath Linear Regression

kmath_linear_regression_hc.kt

import org.ojalgo.random.Normal
import scientifik.kmath.linear.MatrixContext
import scientifik.kmath.linear.VirtualMatrix
import scientifik.kmath.structures.Matrix
import java.net.URL
import kotlin.math.pow

fun main() {

    val points = URL("https://tinyurl.com/y58sesrr")
            .readText().split(Regex("\\r?\\n"))
            .asSequence()
            .drop(1)
            .filter { it.isNotEmpty() }
            .map { it.split(",").map { it.toDouble() }.toDoubleArray() }
            .toMatrix()

    val x = points.extractColumn(0)
    val y = points.extractColumn(1)

    val normal = Normal(0.0,1.0)

    val epochs = 200_000

    val n = x.rowNum.toDouble()

    var m = 0.0
    var b = 0.0

    var bestLoss = 10000000000000.0

    repeat(epochs){
        val mAdjust = normal.get()
        val bAdjust = normal.get()

        m += mAdjust
        b += bAdjust

        val newLoss = ((y - (x * m + b)).square() ).sum() / n

        if (newLoss < bestLoss) {
            bestLoss = newLoss
        } else {
            m -= mAdjust
            b -= bAdjust
        }
    }

    println("f(x) = ${m}x + $b")
}


fun Iterable<DoubleArray>.toMatrix() = toList().let {
    MatrixContext.real.produce(it.size,it[0].size) {row,col -> it[row][col] }
}

fun Sequence<DoubleArray>.toMatrix() = toList().let {
    MatrixContext.real.produce(it.size,it[0].size) {row,col -> it[row][col] }
}

fun <T,N: Any> Iterable<T>.toVirtualMatrix(vararg fieldMappers: (T) -> N) = toList().let { backingList ->
    VirtualMatrix(rowNum = backingList.size, colNum = fieldMappers.count()) { row,col ->
        fieldMappers[col](backingList[row])
    }
}

fun <T,N: Any> Sequence<T>.toVirtualMatrix(vararg fieldMappers: (T) -> N) = toList().toVirtualMatrix(*fieldMappers)

operator fun Matrix<Double>.times(double: Double) = VirtualMatrix(rowNum,colNum) { row,col ->
    this@times[row,col] * double
}

fun Matrix<Double>.square() =  MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@square[row,col].let { it.pow(2) }
}

operator fun Matrix<Double>.plus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@plus[row,col] + double
}

operator fun Matrix<Double>.minus(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@minus[row,col] - double
}

operator fun Matrix<Double>.div(double: Double) = MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@div[row,col] / double
}

operator fun Double.times(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
    matrix[row,col] * this
}

operator fun Double.plus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
    matrix[row,col] + this
}

operator fun Double.minus(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
    matrix[row,col] - this
}

operator fun Double.div(matrix: Matrix<Double>) = MatrixContext.real.produce(matrix.rowNum, matrix.colNum) { row, col ->
    matrix[row,col] / this
}


operator fun Matrix<Double>.times(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@times[row,col] * other[row,col]
}

operator fun Matrix<Double>.plus(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@plus[row,col] + other[row,col]
}

operator fun Matrix<Double>.minus(other: Matrix<Double>) = MatrixContext.real.produce(rowNum, colNum) { row, col ->
    this@minus[row,col] - other[row,col]
}

fun Matrix<Double>.extractColumn(columnIndex: Int) = extractColumns(columnIndex..columnIndex)

fun Matrix<Double>.extractColumns(columnRange: IntRange) = MatrixContext.real.produce(rowNum, columnRange.count()) { row, col ->
    this@extractColumns[row, columnRange.start + col]
}

fun Matrix<Double>.extractRow(rowIndex: Int) = extractRows(rowIndex..rowIndex)

fun Matrix<Double>.extractRows(rowRange: IntRange) = MatrixContext.real.produce(rowRange.count(), colNum) { row, col ->
    this@extractRows[rowRange.start + row, col]
}

fun Matrix<Double>.appendRight(other: Matrix<Double>) =  MatrixContext.real.produce(rowNum + other.rowNum, colNum) { row, col ->
    when {
        col < this@appendRight.colNum -> this@appendRight[row,col]
        else -> other[row, col - this@appendRight.colNum]
    }
}

fun Matrix<Double>.appendBottom(other: Matrix<Double>) = MatrixContext.real.produce(rowNum + other.rowNum, colNum) { row, col ->

    when {
        row < this@appendBottom.rowNum -> this@appendBottom[row,col]

        else -> other[row - this@appendBottom.rowNum, col]
    }
}

fun Matrix<Double>.sum() = this.elements().map { (dim,value) -> value  }.sum()