Почему MATLAB так быстро умножает матрицы?

этот вопрос повторяется и должен быть дан более четкий ответ, чем "Matlab использует высокооптимизированные библиотеки" или "Matlab использует MKL" на один раз на Stackoverflow.

история:

матричное умножение (вместе с матрично-векторным, векторно-векторным умножением и многими матричными разложениями) является (являются) наиболее важными проблемами в линейной альгребре. Инженеры решают эти проблемы с компьютерами с самого начала дни.

Я не эксперт по истории, но, видимо, тогда все просто переписали свою версию Fortran с помощью простых циклов. Затем появилась некоторая стандартизация с идентификацией "ядер" (основных подпрограмм), которые необходимы для решения большинства задач линейной алгебры. Эти основные операции были затем стандартизированы в спецификации под названием: Основные подпрограммы линейной алгебры (BLAS). Затем инженеры могли бы назвать эти стандартные, хорошо протестированные процедуры BLAS в своем коде, что делает их работу намного проще.

Блас:

BLAS эволюционировал от уровня 1 (первая версия, которая определяла скалярно-векторные и векторно-векторные операции) до уровня 2 (векторно-матричные операции) до уровня 3 (матрично-матричные операции) и обеспечивал все больше и больше "ядер", поэтому стандартизировал все больше и больше фундаментальных операций линейной алгебры. Оригинальные реализации Fortran 77 по-прежнему доступны на библиотека netlib по сайт.

к лучшей производительности:

таким образом, на протяжении многих лет (особенно между выпусками BLAS уровня 1 и уровня 2: начало 80-х годов) аппаратное обеспечение изменилось с появлением векторных операций и иерархий кэша. Эти изменения позволили существенно повысить производительность подпрограмм BLAS. Различные поставщики тогда пришли вместе с их реализацией процедур BLAS, которые были все более и более эффективными.

I не знаю всех исторических реализаций (я тогда не родился и не был ребенком), но два из самых заметных вышли в начале 2000-х годов: Intel MKL и GotoBLAS. Ваш Matlab использует Intel MKL, который является очень хорошим, оптимизированным BLAS, и это объясняет отличную производительность, которую вы видите.

технические детали на умножении матрицы:

так почему же Matlab (MKL) так быстро в dgemm (общая матрица двойной точности-матричное умножение)? В простыми словами: потому что он использует векторизацию и хорошее кэширование данных. В более сложных терминах: см. статьи предоставлено Джонатаном Муром.

в принципе, когда вы выполняете свое умножение в коде C++, который вы предоставили, вы совсем не дружелюбны к кэшу. Поскольку я подозреваю, что вы создали массив указателей на массивы строк, ваш доступ во внутреннем цикле к k-му столбцу "matice2":matice2[m][k] очень медленно. Действительно, когда вы получаете доступ matice2[0][k], вы должны получить K-й элемент массива 0 вашей матрицы. Затем в следующей итерации вы должны получить доступ к matice2[1][k], который является k-м элементом другого массива (массив 1). Затем на следующей итерации вы получаете доступ к еще одному массиву и так далее... Так как вся матрица matice2 не может поместиться в самые высокие кэши (это 8*1024*1024 bytes large), программа должна извлечь нужный элемент из основной памяти, теряя много времени.

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

timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
    for (int q = 0; q < rozmer; q++)
    {
        tempmat[p][q] = matice2[q][p];
    }
}
for(int j = 0; j < rozmer; j++)
{
    for (int k = 0; k < rozmer; k++)
    {
        temp = 0;
        for (int m = 0; m < rozmer; m++)
        {
            temp = temp + matice1[j][m] * tempmat[k][m];
        }
        matice3[j][k] = temp;
    }
}
timer.stop();

таким образом, вы можете видеть, как просто локальность кэша значительно увеличила производительность вашего кода. Теперь реально dgemm реализации используют это на очень обширном уровне: они выполняют умножение на блоках матрицы, определяемой размером TLB (Translation lookaside buffer, long story short: what может эффективно кэшироваться), так что они передают процессору ровно тот объем данных, который он может обрабатывать. Другой аспект-векторизация, они используют векторизованные инструкции процессора для оптимальной пропускной способности команд, что вы действительно не можете сделать из своего кросс-платформенного кода C++.

наконец, люди, утверждающие, что это из–за алгоритма Штрассена или Копперсмита-виноград, ошибаются, оба эти алгоритма не реализуемы на практике из-за аппаратного обеспечения соображения, упомянутые выше.

вот почему. MATLAB не выполняет наивное матричное умножение, перебирая каждый элемент так, как вы делали в своем коде на C++.

конечно, я предполагаю, что вы просто использовать C=A*B вместо того, чтобы писать функцию умножения самостоятельно.

Matlab включил LAPACK некоторое время назад, поэтому я предполагаю, что их матричное умножение использует что-то, по крайней мере, так быстро. Исходный код и документация LAPACK легко доступны.

вы также можете посмотреть на статью Гото и Ван де Гейна " Анатомия высокопроизводительной матрицы Умножение " на http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.140.1785&rep=rep1&type=pdf

при умножении матрицы вы используете наивный метод умножения, который занимает время O(n^3).

существует алгоритм умножения матриц, который принимает O(n^2.4). Что означает, что в n=2000 ваш алгоритм требует ~100 раз больше вычислений, чем лучший алгоритм.
Вы действительно должны проверить страницу Википедии для умножения матрицы для получения дополнительной информации об эффективных способах ее реализации.

вы должны быть осторожны, делая справедливые сравнения с C++. Можете ли вы опубликовать код C++, который показывает основные внутренние циклы, которые вы используете для умножения матрицы? В основном, меня беспокоит ваш макет памяти и делаете ли вы что-то расточительно.

Я написал c++ матричное умножение, которое так же быстро, как и Matlab, но оно было немного осторожным. (EDIT: до того, как Matlab использовал для этого графические процессоры.)

вы можете быть практически уверены, что Matlab тратит впустую очень мало циклов на этих "встроенных" функциях. Мой вопрос в том, где вы тратите циклов? (Без обид)

ответ LAPACK и Блас библиотеки делают MATLAB ослепительно быстрым в матричных операциях, а не в каком-либо собственном коде людей в MATLAB.

использовать LAPACK и/или Блас библиотеки в вашем коде C++ для матричных операций, и вы должны получить аналогичную производительность, как MATLAB. Эти библиотеки должны быть в свободном доступе на любой современной системе и части были разработаны в течение десятилетий в академических кругах. Обратите внимание, что есть несколько реализаций, включая некоторые закрытые исходники, такие как Intel MKL.

обсуждение того, как BLAS получает высокую производительность здесь.

---

кстати, это серьезная боль в моем опыте, чтобы вызвать библиотеки LAPACK непосредственно из c (но стоит того). Вам нужно очень точно прочитать документацию.

в зависимости от вашей версии Matlab, я считаю, что он уже может использовать ваш GPU.

другое дело; Matlab отслеживает многие свойства вашей матрицы; имеет свою диагональ, герметичность и т. д. и специализируется на своих алгоритмах, основанных на них. Может быть, его специализация основана на нулевой матрице, которую вы передаете, или что-то в этом роде? Может быть, это кэширование повторяющихся вызовов функций, что портит ваши тайминги? Возможно, он оптимизирует повторяющуюся неиспользуемую матрицу продукты?

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

использование двойников и одного сплошного массива вместо трех отдельных приводит мой код C# к почти тем же результатам, что и C++/Java (с вашим кодом: 1024 - было немного быстрее, 2048 - было около 140s и 4096 - было около 22 минут)

                1024x1024   2048x2048   4096x4096
                ---------   ---------   ---------
your C++ (ms)   6137.10     64369.29     551390.93
my C# (ms)      9730.00     90875.00    1062156.00

вот мой код:

    const int rozmer = 1024;
    double[][] matice1 = new double[rozmer * 3][];
    Random rnd = new Random();

    public Form1()
    {
        InitializeComponent();

        System.Threading.Thread thr = new System.Threading.Thread(new System.Threading.ThreadStart(() =>
        {

            string res = "";
            Stopwatch timer = new Stopwatch();
            timer.Start();

            double temp = 0;
            int r2 = rozmer * 2;

            for (int i = 0; i < rozmer*3; i++)
            {
                if (matice1[i] == null)
                {
                    matice1[i] = new double[rozmer];
                    {
                        for (int e = 0; e < rozmer; e++)
                        {
                            matice1[i][e] = rnd.NextDouble();
                        }
                    }
                }
            }
            timer.Stop();
            res += timer.ElapsedMilliseconds.ToString();

            int j = 0; int k = 0; int m = 0;

            timer.Reset();
            timer.Start();
            for (j = 0; j < rozmer; j++)
            {
                for (k = 0; k < rozmer; k++)
                {
                    temp = 0;
                    for (m = 0; m < rozmer; m++)
                    {
                        temp = temp + matice1[j][m] * matice1[m + rozmer][k];
                    }
                    matice1[j + r2][k] = temp;
                }
            }
            timer.Stop();
            this.Invoke((Action)delegate
            {
                this.Text = res + " : " + timer.ElapsedMilliseconds.ToString();
            });
        }));
        thr.Start();
    }

вы проверили, что все реализации использовали многопоточные оптимизации для алгоритма ? И они использовали один и тот же алгоритм умножения ?

Я действительно сомневаюсь в этом.

Matlab не является по своей сути быстрым, вы, вероятно, использовали медленные реализации.

алгоритмы эффективного умножения матриц

общий ответ на вопрос" почему matlab быстрее делает xxx, чем другие программы " заключается в том, что matlab имеет много встроенных, оптимизированных функций.

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

Это можно интерпретировать двумя способами:

1) общий / теоретический способ: Matlab не значительно быстрее, вы просто делаете эталон неправильный

2) реалистичный способ: для этого материала Matlab быстрее на практике, потому что языки, такие как c++, Слишком легко используются неэффективными способами.

MATLAB использует высоко оптимизированную реализацию LAPACK от Intel, известную как Библиотека Ядра Intel Math (Intel MKL) - конкретно dgemm функция. Скорость эта библиотека использует в своих интересах функции процессора, включая инструкции SIMD и многоядерные процессоры. Они не документируют, какой конкретный алгоритм они используют. Если вы должны были вызвать Intel MKL из C++ , вы должны увидеть аналогичную производительность.

Я не уверен, что библиотека MATLAB использует для Умножение GPU, но, вероятно, что-то вроде nVidia CUBLAS.

это медленно в C++, потому что вы не используете многопоточность. По существу, если A = B C, где все они являются матрицами, первая строка A может быть вычислена независимо от второй строки и т. д. Если A, B и C все N на N матриц, вы можете ускорить умножение на коэффициент N^2, как

a_{i, j} = sum_{k} b_{i,k} c_{k, j}

Если вы используете, скажем, Eigen [ http://eigen.tuxfamily.org/dox/GettingStarted.html ], многопоточность встроена и количество потоков регулируется.

резкий контраст обусловлен не только удивительной оптимизацией Matlab (как уже обсуждалось многими другими ответами), но и тем, как вы сформулировали матрицу как объект.

похоже, вы сделали матрицу списком списков? Список списков содержит указатели на списки, которые затем содержат элементы матрицы. Расположение содержащихся в них списков назначается произвольно. Как вы перебираете свой первый индекс (номер строки?), время доступа к памяти очень значительно. В сравнение, почему бы вам не попробовать реализовать матрицу в виде одного списка/вектора, используя следующий метод?

#include <vector>

struct matrix {
    matrix(int x, int y) : n_row(x), n_col(y), M(x * y) {}
    int n_row;
    int n_col;
    std::vector<double> M;
    double &operator()(int i, int j);
};

и

double &matrix::operator()(int i, int j) {
    return M[n_col * i + j];
}

тот же алгоритм умножения должен использоваться так, чтобы число флопов было одинаковым. (н^3 для квадратных матриц размера n)

Я прошу вас рассчитать время так, чтобы результат был сопоставим с тем, что вы имели ранее (на той же машине). С помощью сравнения вы точно покажете, насколько значительным может быть время доступа к памяти быть!