Регистрация

[JavaScript, Программирование, Клиентская оптимизация, Математика] Кэширование данных увеличивает скорость даже в неожиданных случаях

Ответить на тему
Автор Сообщение
news_bot ®

Стаж: 6 лет 9 месяцев
Сообщений: 27286

Создавать темы news_bot ® написал(а)
13-Апр-2021 18:31
Цитировать

Нас учат, что чтение данных из оперативной памяти — ужасно долгая операция. Приводят аналогии с офисом и удалённым складом, заставляют писать 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().Помимо wholeMod в вычитании строк матриц, используется . Кроме того, появилась вышеуказанная проблема с ограничением mulValue. Во всех этих случаях деление описывается формулой 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.
 ===========
 Источник:
habr.com
===========
Похожие новости: Теги для поиска: #_javascript, #_programmirovanie (Программирование), #_klientskaja_optimizatsija (Клиентская оптимизация), #_matematika (Математика), #_optimizatsija (оптимизация), #_javascript, #_matritsy (матрицы), #_proizvoditelnost (производительность), #_udivitelnoe (удивительное), #_javascript, #_programmirovanie (
Программирование
), #_klientskaja_optimizatsija (
Клиентская оптимизация
), #_matematika (
Математика
)
Профиль  ЛС 
Показать сообщения:     
Ответить на тему

Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы

Текущее время: 22-Ноя 13:41
Часовой пояс: UTC + 5