Библиотека работы с вещественными матрицами

Опись имущества:

1. class D2DoubleMatrix: содержит представление двумерной вещественной матрицы

   1. D2DoubleMatrix.getEntryArray() возвращает DoubleArray(), в котором содержатся все элементы матрицы, вытянутые в линию строчка за строчкой

   2. D2DoubleMatrix.shape() возвращает пару Pair<Int, Int>, где Pair.first количество строк матрицы, Pair.second количество столбцов

   3. D2DoubleMatrix.numberAt(r: Int, c: Int) элемент в строке r и столбце c

   4. D2DoubleMatrix.setAt(r: Int, c: Int, val: Double) установка значения в строке r столбце c

   5. D2DoubleMatrix.toString() красивый вывод прилагается! (совсем красивый вывод ещё надо доработать, чтобы точки были под дочками, запятые под запятыми)

   6. == сравнивает матрицы на равенство с точностью до 1e-10 (может быть не самое лучшее решение)

   7. *, -, +. Также есть умножение на скаляр справа.

2. abstract class MatrixCreator содержит фабричные методы создния матриц

   1. fromRows(vararg rows: DoubleArray) создает матицу из последовательности строк

   2. zero(row: Int, column: Int) матрица, заполненная нулями, в которой row строк column столбцов

   3. randomElementUniformMatrix(row: Int, column: Int, leftBound: Double, rightBound: Double) матрица формы (row, column), каждый элемент которой случайня величина равонмерно распределенная на [leftBound, rightBound]

   4. clone(A: D2DoubleMatrix) возвращает копию матрицы

   5. transpose(A: D2DoubleMatrix) возвращает транспонированную матрицу

   6. fromFlatArray(row: Int, column: Int, arr: DoubleArray) создает матрицу формы (row, column)из плоского массива. Должно быть выполнено arr.size() == row * column

   7. embed(A: D2DoubleMatrix, rows: Int, columns: Int) вкладывает матрицу в левый верхний угол большей нулевой матрицы.

   8. fun slice(A: D2DoubleMatrix, fromRowInclusive: Int, toRowExclusive: Int, fromColumnInclusive: Int, toColumnExclusive: Int ) вырезает из матрицы подматрицу

3. interface multiplyStrategy общий класс для всех алгоритмов умножния матрицы на матрицу. Далее перечислены реализации:

   1. DummyMultiplier обычное умножение двух матриц по определению

   2. DummyCacheFriendly то же что и обычное, но старается бегать по подряд идущим элементам массивов.

   3. Strassen реализация алгоритма Штрассена. Параметры внутри класса, которые можно было бы настраивать, если бы я не объявил их приватными:

      1. private strassenThreshold = 50. Если хоть одна размерность перемножаемых матриц становиться меньше этого значения, алгоритм переключается на defaultMultiplier

      2. private defaultMultiplier какой-то алгоритм, асимптотически худший, чем алгоритм Штрассена, но в котором константа лучше

   4. SmallMatrixMultiplier класс, чтобы умножать маленькие матрицы. Сам выбирает, что надо использовать: DummyMultiplier или DummyCacheFriendly

   Сравнительная таблица скорости различных multiplyStrategy (тесты проведены на моей локальной машине):

   Mat(m, n) x Mat(n, k)	DummyMultiplier	DummyCacheFriendly	Strassen
   (2026, 2026)x(2026, 2026)	53 s	14 s	13 s
   (100, 100000)x(100000, 100)	40 s	20 s	20 s

4. user.kt примеры кода

5. в разделе src/test можно найти некоторые (достаточно скудные) юнит тесты

Примеры использования

fun main() {
    val mt1 = MatrixCreator.fromRows(doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0))
    val mt2 = MatrixCreator.transpose(mt1)
    println(mt1)
    /*
    *[[1.0, 2.0, 3.0, 4.0, 5.0, 6.0]]
    */
    println(mt2)
    /*
    * [[1.0],
    *  [2.0],
    *  [3.0],
    *  [4.0],
    *  [5.0],
    *  [6.0]]
    */
    println(mt1 * mt2)
    /*
    * [[91.0]]
    */
    println(mt2 * mt1)
    /*
    * [[1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
    *  [2.0, 4.0, 6.0, 8.0, 10.0, 12.0],
    *  [3.0, 6.0, 9.0, 12.0, 15.0, 18.0],
    *  [4.0, 8.0, 12.0, 16.0, 20.0, 24.0],
    *  [5.0, 10.0, 15.0, 20.0, 25.0, 30.0],
    *  [6.0, 12.0, 18.0, 24.0, 30.0, 36.0]]
    */
    println(mt2 * mt1 * mt2)
    /*
    * [[91.0],
    * [182.0],
    * [273.0],
    * [364.0],
    * [455.0],
    * [546.0]]
    */

    println("we are going to multiply something big")

    //большие случайные матрицы
    val mt3 = MatrixCreator.randomElementUniformMatrix(1000, 1000, 0.0, 1.0)
    val mt4 = MatrixCreator.randomElementUniformMatrix(1000, 1000, 0.0, 1.0)

    val mt5: D2DoubleMatrix
    val mt6: D2DoubleMatrix

    var timeDelta = kotlin.system.measureTimeMillis {
        mt5 = mt3 * mt4
    }
    println("time elapsed: ${timeDelta / 1000.0} s")
    // time elapsed: 2.837 s
    
    println("now do it manually")
    timeDelta = kotlin.system.measureTimeMillis {
        mt6 = MatrixCreator.zero(1000, 1000)
        for (i in 0 until mt6.shape().first) {
            for (j in 0 until mt6.shape().second) {
                var aij = 0.0
                for (k in 0 until mt4.shape().first) {
                    aij += mt3.numberAt(i, k) * mt4.numberAt(k, j)
                }
                mt6.setAt(i, j, aij)
            }

        }
    }
    println("time elapsed: ${timeDelta / 1000.0}")
    // time elapsed: 10.49
    
    println("matrix is too large, but we can pick little piece")

    println(MatrixCreator.slice(mt5, 100, 110, 100, 110))
    println(MatrixCreator.slice(mt6, 100, 110, 100, 110))
    /*
    * [[246.060585504873, 250.16693175308245, 240.68578746445698, 255.29232309450074, 253.9489627454569],
    *  [254.45035890564046, 247.31633774856314, 248.30619847906286, 251.8215115935202, 254.38594309588066],
    *  [246.3234365478975, 245.8108942512048, 245.469841326687, 244.97957579338515, 249.9707712454963],
    *  [246.28611172880207, 247.5417773235133, 241.52132635472313, 248.00274049050913, 249.73343167845655],
    *  [257.1752544076303, 255.68853949321317, 249.18423743491036, 255.2947615556086, 254.46422883341526]]
    */
    // тут я хотел считерить и скопировать то что выше, но потом заметил
    // что результаты отличаются в 11 знаке
    /*
    * [[246.06058550487347, 250.16693175308234, 240.68578746445627, 255.2923230945008, 253.94896274545616],
    *  [254.45035890563992, 247.3163377485639, 248.30619847906314, 251.82151159351898, 254.3859430958808],
    *  [246.3234365478975, 245.81089425120524, 245.46984132668643, 244.97957579338456, 249.97077124549608],
    *  [246.28611172880161, 247.54177732351323, 241.5213263547229, 248.0027404905089, 249.7334316784566],
    *  [257.17525440763046, 255.6885394932124, 249.18423743490993, 255.2947615556086, 254.46422883341518]]
    */
}

Послесловие

1. Конечно, это должна быть библиотека умножения матриц любого наследника Number, а не умножения вещественных числел. Код переносится на натуральные с точностью до замены по всему проекту слова Double на Int. Но я быстро не смог придумать как реализовать. (Менее убедительное оправдание: а вот в js вообще все числа с плавающей точкой, а я чем хуже)

2. Надо дальше оптимизировать алгоритм Штрассена, так как пока он совсем не впечатляет. Или на таких маленьких матрицах (которые влезают в RAM) его сложно ощутить.

3. Нужно добавить умножение разреженных матриц.

4. Пока что всё однопоточное.

5. Если v это очень длинная строка из нулей, то v * transpose(v) * v зовет переполнение памяти из-за левой ассоциативности, когда как если расставить скобки правильно, то ничего подобоного не возникнет. Это тоже надо бы учесть