13.04.2021 12:49
Sixshaman 17 4000 Источник

Нас учат, что чтение данных из оперативной памяти — ужасно долгая операция. Приводят аналогии с офисом и удалённым складом, заставляют писать cache-friendly код и внушают смертельный страх перед промахами кэша. Ещё нас учат, что процессоры отлично умеют считать числа, и часто быстрее вычислить результат дважды, чем сохранять его в памяти. Оказывается, это не всегда так.

Эта статья основана на реальном проекте и реальном коде, который был ускорен с помощью кэша почти в полтора раза. Весь код написан на JavaScript.

Задача

Допустим, у нас есть матрица A порядка 2000x2000. Нужно посчитать обратную ей матрицу по простому модулю N. Другими словами, надо найти такую матрицу A-1, что AA-1 mod N = E.

Поскольку вычисления у нас происходят в поле по модулю, итерационные методы обращения нам не подойдут. Будем использовать старый добрый метод Гаусса.

Этот пост посвящён оптимизации метода Гаусса под данный конкретный случай. В реальном проекте вычисление обратной матрицы происходит в отдельном WebWorker, в данном примере обойдёмся главным потоком.

Вспомогательные функции

Для работы программы нам потребуется четыре вспомогательные функции. Первая — вычисление (1 / x) mod N по расширенному алгоритму Евклида:

Код
function invModGcdEx(x, domain)
{
    if(x === 1)
    {
        return 1;
    }
    else
    {
        //В случае 0 или делителя нуля возвращается 0, означающий некий "некорректный результат"
        if(x === 0 || domain % x === 0)
        {
            return 0;
        }
        else
        {
            //Расширенный алгоритм Евклида, вычисляющий такое число tCurr, что tCurr * x + rCurr * N = 1
            //Другими словами, существует такое число rCurr, при котором tCurr * x mod N = 1
            let tCurr = 0;
            let rCurr = domain;
            let tNext = 1;
            let rNext = x;
            while(rNext !== 0)
            {
                let quotR = Math.floor(rCurr / rNext);
                let tPrev = tCurr;
                let rPrev = rCurr;

                tCurr = tNext;
                rCurr = rNext;

                tNext = Math.floor(tPrev - quotR * tCurr);
                rNext = Math.floor(rPrev - quotR * rCurr);
            }

            tCurr = (tCurr + domain) % domain;
            return tCurr;
        }
    }
}

Вторая — корректное целочисленное деление по модулю. Наивное вычисление c = a % b во всех языках программирования не будет давать математически верный результат, если a — отрицательное число. Поэтому заведём функцию, которая будет делить правильно:

Код
function wholeMod(x, domain)
{
    return ((x % domain) + domain) % domain;
}

Последние две функции относятся к операциям над строками матрицы. Первая — вычитание из строки матрицы домноженной на число другой:

Код
function mulSubRow(rowLeft, rowRight, mulValue, domain)
{
    for(let i = 0; i < rowLeft.length; i++)
    {
        rowLeft[i] = wholeMod(rowLeft[i] - mulValue * rowRight[i], domain);
    }
}

Последняя нужная нам функция — умножение строки матрицы на число:

Код
function mulRow(row, mulValue, domain)
{
    for(let i = 0; i < row.length; i++)
    {
        row[i] = (row[i] * mulValue) % domain;
    }
}

Обращение матрицы

Начнём с обычной наивной реализации. Создаём единичную матрицу, проходим прямым ходом, потом проходим обратным. На каждом шаге производим одинаковые операции над исходной матрицей и над только что созданной единичной.

Код
function invertMatrix(matrix, domain)
{
    let matrixSize = matrix.length;

    //Инициализируем обратную матрицу единичной
    let invMatrix = [];
    for(let i = 0; i < matrixSize; i++)
    {
        let matrixRow = new Uint8Array(matrixSize);
        matrixRow.fill(0);

        matrixRow[i] = 1;
        invMatrix.push(matrixRow);
    }

    //Прямой ход: приведение матрицы к ступенчатому виду
    for(let i = 0; i < matrixSize; i++)
    {
        let thisRowFirst = matrix[i][i];
        if(thisRowFirst === 0 || (thisRowFirst !== 1 && domain % thisRowFirst === 0)) //Первый элемент строки 0 или делитель нуля, меняем строку местами со следующей строкой, у которой первый элемент не 0
        {
            for(let j = i + 1; j < matrixSize; j++)
            {
                let otherRowFirst = matrix[j][i];
                if(otherRowFirst !== 0 && (otherRowFirst === 1 || domain % otherRowFirst !== 0)) //Нашли строку с ненулевым первым элементом
                {
                    thisRowFirst = otherRowFirst;
                    
                    let tmpMatrixRow = matrix[i];
                    matrix[i]        = matrix[j];
                    matrix[j]        = tmpMatrixRow;

                    let tmpInvMatrixRow = invMatrix[i];
                    invMatrix[i]        = invMatrix[j];
                    invMatrix[j]        = tmpInvMatrixRow;

                    break;
                }
            }
        }

        //Обнуляем первые элементы всех строк после первой, отнимая от них (otherRowFirst / thisRowFirst) * x mod N
        let invThisRowFirst = invModGcdEx(thisRowFirst, domain);
        for(let j = i + 1; j < matrixSize; j++)
        {
            let otherRowFirst = matrix[j][i];
            let mulValue      = invThisRowFirst * otherRowFirst;

            if(otherRowFirst !== 0 && (otherRowFirst === 1 || domain % otherRowFirst !== 0))
            {
                mulSubRow(matrix[j],    matrix[i],    mulValue, domain);
                mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
            }
        }
    }

    //Обратный ход - обнуление всех элементов выше главной диагонали
    let matrixRank = matrixSize;
    for(let i = matrixSize - 1; i >= 0; i--)
    {
        let thisRowLast    = matrix[i][i];
        let invThisRowLast = invModGcdEx(thisRowLast, domain);
        for(let j = i - 1; j >= 0; j--)
        {
            let otherRowLast = matrix[j][i];
            let mulValue     = invThisRowLast * otherRowLast;

            if(otherRowLast !== 0 && (otherRowLast === 1 || domain % otherRowLast !== 0))
            {
                mulSubRow(matrix[j],    matrix[i],    mulValue, domain);
                mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
            }
        }

        if(thisRowLast !== 0 && domain % thisRowLast !== 0)
        {
            mulRow(matrix[i],    invThisRowLast, domain);
            mulRow(invMatrix[i], invThisRowLast, domain);
        }

        if(matrix[i].every(val => val === 0))
        {
            matrixRank -= 1;
        }
    }

    return {inverse: invMatrix, rank: matrixRank};
}

Проверим скорость на матрице 500 x 500, заполненной случайными значениями из поля Z / 29. После 5 испытаний получаем среднее время выполнения в ~9.4с. Можем ли мы сделать лучше?

Первое, что мы можем заметить — в поле Z / N не больше N обратных элементов. Чтобы избежать многократного вызова алгоритма Евклида, мы можем вычислить все обратные значения один раз и при надобности брать уже готовые. Изменим функцию соответствующим образом:

Код
function invertMatrixCachedInverses(matrix, domain)
{
    let matrixSize = matrix.length;

    //Инициализируем обратную матрицу единичной
    let invMatrix = [];
    for(let i = 0; i < matrixSize; i++)
    {
        let matrixRow = new Uint8Array(matrixSize);
        matrixRow.fill(0);

        matrixRow[i] = 1;
        invMatrix.push(matrixRow);
    }

    //Вычисляем все обратные элементы заранее
    let domainInvs = [];
    for(let d = 0; d < domain; d++)
    {
        domainInvs.push(invModGcdEx(d, domain));
    }

    //Прямой ход: приведение матрицы к ступенчатому виду
    for(let i = 0; i < matrixSize; i++)
    {
        let thisRowFirst = matrix[i][i];
        if(domainInvs[thisRowFirst] === 0) // <--- Первый элемент строки 0 или делитель нуля, меняем строку местами со следующей строкой, у которой первый элемент не 0
        {
            for(let j = i + 1; j < matrixSize; j++)
            {
                let otherRowFirst = matrix[j][i];
                if(domainInvs[otherRowFirst] !== 0) // <--- Нашли строку с ненулевым первым элементом
                {
                    thisRowFirst = otherRowFirst;
                    
                    let tmpMatrixRow = matrix[i];
                    matrix[i]        = matrix[j];
                    matrix[j]        = tmpMatrixRow;

                    let tmpInvMatrixRow = invMatrix[i];
                    invMatrix[i]        = invMatrix[j];
                    invMatrix[j]        = tmpInvMatrixRow;

                    break;
                }
            }
        }

        //Обнуляем первые элементы всех строк после первой, отнимая от них (otherRowFirst / thisRowFirst) * x mod N
        let invThisRowFirst = domainInvs[thisRowFirst]; // <---
        for(let j = i + 1; j < matrixSize; j++)
        {
            let otherRowFirst = matrix[j][i];
            let mulValue      = invThisRowFirst * otherRowFirst;

            if(domainInvs[otherRowFirst] !== 0) // <---
            {
                mulSubRow(matrix[j],    matrix[i],    mulValue, domain);
                mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
            }
        }
    }

    //Обратный ход - обнуление всех элементов выше главной диагонали
    let matrixRank = matrixSize;
    for(let i = matrixSize - 1; i >= 0; i--)
    {
        let thisRowLast    = matrix[i][i];
        let invThisRowLast = domainInvs[thisRowLast]; // <---
        for(let j = i - 1; j >= 0; j--)
        {
            let otherRowLast = matrix[j][i];
            let mulValue     = invThisRowLast * otherRowLast;

            if(domainInvs[otherRowLast] !== 0) // <---
            {
                mulSubRow(matrix[j],    matrix[i],    mulValue, domain);
                mulSubRow(invMatrix[j], invMatrix[i], mulValue, domain);
            }
        }

        if(domainInvs[thisRowLast] !== 0) // <---
        {
            mulRow(matrix[i],    invThisRowLast, domain);
            mulRow(invMatrix[i], invThisRowLast, domain);
        }

        if(matrix[i].every(val => val === 0))
        {
            matrixRank -= 1;
        }
    }

    return {inverse: invMatrix, rank: matrixRank};
}

Замерим на тех же условиях и получаем результат в те же ~9.4с. Прироста нет, потому что даже при относительно долгом вычислении алгоритма Евклида он вычисляется всего один раз для каждой строки матрицы и особого вклада во время не приносит. Замерим производительность и посмотрим, что ещё можно улучшить.

72% времени занимает деление по модулю при сложении строк матрицы! Ну что тут сказать, деление по модулю, пусть и немного модифицированное для отрицательных чисел — это элементарная операция и ускорять её некуда. Алгоритм поменять тоже не получится, из чего мы делаем вывод, что дальнейшее улучшение невозможно и статью можно закрывать.

...Или всё же возможно?

Если деление по модулю занимает столько времени, может, все возможные результаты тоже стоит сохранить в кэш? Даже если это не поможет, попытаться всё равно стоит — при текущем времени выполнения функция неюзабельна.

Итак, используется wholeMod()только в функции mulSubRow():

rowLeft[i] = wholeMod(rowLeft[i] - mulValue * rowRight[i], domain);

Нам нужно для всех возможных значений x = a - b * c в поле Z / N сохранить результат выражения x mod N. Воспользоваться периодичностью мы не сможем, потому что тогда для вычисления индекса снова придётся использовать деление по модулю. В итоге при 0 <= a, b, c < N получаем N + (N - 1)^2 возможных значений. Много, но деваться некуда.

Из этих значений (N - 1)^2 значений меньше 0. Поскольку отрицательные индексы невозможны, при индексировании значением a - b * c к нему нужно прибавить (N - 1)^2. Тогда функция для сложения строк модифицируется:

Код
function mulSubRowCached(rowLeft, rowRight, mulValue, wholeModCache, cacheIndexOffset)
{
    for(let i = 0; i < rowLeft.length; i++)
    {
        rowLeft[i] = wholeModCache[rowLeft[i] - mulValue * rowRight[i] + cacheIndexOffset];
    }
}

Заметим, что эта функция накладывает ограничение на mulValue — его значение не может быть больше domain и перед вызовом функции его тоже надо привести в наше поле Z / N. Кроме этого, обычное деление по модулю используется в функции mulRow(). Во всех этих случаях деление описывается формулой x = (a * b) mod N. Зная, что кэш хранит значения x = (c - a * b) mod N, мы можем вычислить (a * b) mod N, взяв значение кэша при c = 0 и вычтя его из N. Тогда функция для умножения строки на число модифицируется следующим образом:

Код
function mulRowCached(row, mulValue, domain, wholeModCache, cacheIndexOffset)
{
    for(let i = 0; i < row.length; i++)
    {
        row[i] = domain - wholeModCache[cacheIndexOffset - row[i] * mulValue];
    }
}

И получаем новое обращение матрицы:

Код
function invertMatrix(matrix, domain)
{
    let matrixSize = matrix.length;

    //Инициализируем обратную матрицу единичной
    let invMatrix = [];
    for(let i = 0; i < matrixSize; i++)
    {
        let matrixRow = new Uint8Array(matrixSize);
        matrixRow.fill(0);

        matrixRow[i] = 1;
        invMatrix.push(matrixRow);
    }

    //Вычисляем все обратные элементы заранее
    let domainInvs = [];
    for(let d = 0; d < domain; d++)
    {
        domainInvs.push(invModGcdEx(d, domain));
    }

    //Вычисляем кэш деления по модулю
    const сacheIndexOffset = (domain - 1) * (domain - 1);

    let wholeModCache = new Uint8Array((domain - 1) * (domain - 1) + domain); 
    for(let i = 0; i < wholeModCache.length; i++)
    {
        let divisor      = i - сacheIndexOffset;      //[-domainSizeCacheOffset, domainSize - 1]
        wholeModCache[i] = wholeMod(divisor, domain); //Whole mod
    }

    //Прямой ход: приведение матрицы к ступенчатому виду
    for(let i = 0; i < matrixSize; i++)
    {
        let thisRowFirst = matrix[i][i];
        if(domainInvs[thisRowFirst] === 0) //Первый элемент строки 0 или делитель нуля, меняем строку местами со следующей строкой, у которой первый элемент не 0
        {
            for(let j = i + 1; j < matrixSize; j++)
            {
                let otherRowFirst = matrix[j][i];
                if(domainInvs[thisRowFirst] !== 0) //Нашли строку с ненулевым первым элементом
                {
                    thisRowFirst = otherRowFirst;
                        
                    //Меняем строки местами
                    let tmpMatrixRow = matrix[i];
                    matrix[i]        = matrix[j];
                    matrix[j]        = tmpMatrixRow;

                    let tmpInvMatrixRow = invMatrix[i];
                    invMatrix[i]        = invMatrix[j];
                    invMatrix[j]        = tmpInvMatrixRow;

                    break;
                }
            }
        }

        //Обнуляем первые элементы всех строк после первой, отнимая от них (otherRowFirst / thisRowFirst) * x mod N
        let invThisRowFirst = domainInvs[thisRowFirst]; // <---
        for(let j = i + 1; j < matrixSize; j++)
        {
            let otherRowFirst = matrix[j][j];
            if(domainInvs[otherRowFirst] !== 0)
            {
                let mulValue = domain - wholeModCache[сacheIndexOffset - otherRowFirst * invThisRowFirst]; // <---

                mulSubRowCached(matrix[j],    matrix[i],    mulValue, wholeModCache, сacheIndexOffset); // <---
                mulSubRowCached(invMatrix[j], invMatrix[i], mulValue, wholeModCache, сacheIndexOffset); // <---
            }
        }
    }

    //Обратный ход - обнуление всех элементов выше главной диагонали
    let matrixRank = matrixSize;
    for(let i = matrixSize - 1; i >= 0; i--)
    {
        let thisRowLast    = matrix[i][i];
        let invThisRowLast = domainInvs[thisRowLast];

        for(let j = i - 1; j >= 0; j--)
        {
            let otherRowLast = matrix[j][i];
            if(domainInvs[otherRowLast] !== 0)
            {
                let mulValue = domain - wholeModCache[сacheIndexOffset - otherRowLast * invThisRowLast]; // <---

                mulSubRowCached(matrix[j],    matrix[i],    mulValue, wholeModCache, сacheIndexOffset); // <---
                mulSubRowCached(invMatrix[j], invMatrix[i], mulValue, wholeModCache, сacheIndexOffset); // <---
            }
        }

        if(domainInvs[thisRowLast] !== 0)
        {
            mulRowCached(matrix[i],    invThisRowLast, domain, wholeModCache, сacheIndexOffset); // <---
            mulRowCached(invMatrix[i], invThisRowLast, domain, wholeModCache, сacheIndexOffset); // <---
        }

        if(matrix[i].every(val => val === 0))
        {
            matrixRank -= 1;
        }
    }

    return {inverse: invMatrix, rank: matrixRank};
}

Замерим производительность. На той же матрице 500x500 по модулю 29 получаем время выполнения в ~5.4с.

Простите, что?

Нет, серьёзно, как это возможно? Кэшируем результат деления. Операции на два такта. В век супермедленной памяти и супербыстрых процессоров. Получаем прирост в 40%. Как?

Да, использование JavaScript создаёт определённый оверхед. Но JIT его нивелирует. Видимо, либо он нивелирует его недостаточно, либо не всё, чему нас учат про cache-friendly код — правда.

И да, размер кэша растёт квадратично. Но если сравнить среднее время в полях по разному модулю, то прирост будет не сильно отличаться:

В реальном проекте, где был применён этот метод, матрицы не рандомные и прирост ещё заметнее.

Заключение

Можно ли ещё больше ускорить вычисление? К сожалению, больше ни одного способа я не знаю. Я думал над распараллеливанием вычислений, но обращение матриц очень плохо параллелизуется. Поэтому пока остаётся так.

Полный код я выложил на Pastebin.