Сделайте расширение Laplace более эффективным

Question

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

double getDeterminantLaplace(const std::vector<std::vector<double>> vect) {

    int dimension = vect.size();

    if(dimension == 0) {
        return 1;
    }

    if(dimension == 1) {
        return vect[0][0];
    }

    //Formula for 2x2-matrix
    if(dimension == 2) {
        return vect[0][0] * vect[1][1] - vect[0][1] * vect[1][0];
    }

    double result = 0;
    int sign = 1;
    for(int i = 0; i < dimension; i++) {

        //Submatrix
        std::vector<std::vector<double>> subVect(dimension - 1, std::vector<double> (dimension - 1));
        for(int m = 1; m < dimension; m++) {
            int z = 0;
            for(int n = 0; n < dimension; n++) {
                if(n != i) {
                    subVect[m-1][z] = vect[m][n];
                    z++;
                }
            }
        }

        //recursive call
        result = result + sign * vect[0][i] * getDeterminantLaplace(subVect);
        sign = -sign;
    }

    return result;
}

Теперь у меня вопрос: как сделать этот алгоритм более эффективным?

Один из моя идея состоит в том, чтобы не создавать «подматрицы», а просто работать с исходной матрицей, но я действительно не знаю, как это сделать. Что вы думаете об этой идее? Как я могу сделать это в C ++?

У вас есть еще идеи?

Answer 1

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

Следующая оптимизация - это то, что вы уже предложили: не создавать все подматрицы. Вы можете сделать это, создав индексный вектор. Например, если ваша исходная матрица имеет 4 × 4 элемента, вы рекурсируете со следующими индексными векторами:

0: {1, 2, 3}
1: {0, 2, 3}
2: {0, 1, 3}
3: {0, 1, 2}

Вам не нужно каждый раз создавать индексный вектор с нуля: начните с подвектора, который является текущим вектором без его фронта, затем замените i-е место i-й записью текущего убвектора.

При доступе к элементу s[r][c] подматрицы элемент доступа a[r + top][col[c]] из оригинальная матрица. Вы можете определить индекс верхней строки по размерам текущего вектора столбца и исходной матрицы.

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

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

1, 2, 3     0, 2, 3     0, 1, 3     0, 1, 2

Оттуда вы вычислите эти 2 × 2 подматрицы:

   2, 3        2, 3        1, 3        1, 2
1,    3     0,    3     0,    3     0,    2
1, 2        0, 2,       0, 1,       0, 1,  

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

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

В любом случае, вот код, который реализует все оптимизации, которые я описал:

#include <vector>
#include <map>
#include <iostream>

double subdet(const std::vector<std::vector<double> > &a,
    const std::vector<int> &col,
    std::map<std::vector<int>, double> &memo)
{
    int dim = col.size();
    int top = a.size() - dim;

    if (memo.find(col) != memo.end()) {
        return memo[col];
    }   

    if (dim == 2) return a[top + 0][col[0]] * a[top + 1][col[1]]
                       - a[top + 0][col[1]] * a[top + 1][col[0]];

    double result = 0.0;
    int sign = 1;

    std::vector<int> ncol(&col[1], &col[dim]);

    for (int i = 0; i < dim; i++) {
        if (a[top][col[i]]) {
            double d = subdet(a, ncol, memo);

            result = result + sign * a[top][col[i]] * d;
        }

        sign = -sign;

        if (i + 1 < dim) ncol[i] = col[i];
    }

    memo[col] = result;
    return result;
}

double det(const std::vector<std::vector<double> > a)
{
    int dim = a.size();

    if (dim == 0) return 1.0;
    if (dim == 1) return a[0][0];

    std::vector<int> col(dim);
    std::map<std::vector<int>, double> memo;

    for (unsigned i = 0; i < a.size(); i++) col[i] = i;

    return subdet(a, col, memo);
}

Примечания: Карта (двоичное дерево с O (log * 1027) * n ) lookup) действительно должна быть неупорядоченная карта (таблица ha sh с поиском O (1)), но я не смог заставить ее работать, потому что я плохо разбираюсь в C ++. Извините за это.

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

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