Pull to refresh

Снижение сложности вычислений при операциях с векторами и матрицами

Reading time 6 min
Views 7.5K

Введение


Ввиду того, что при решении задач оптимизации, дифференциальных игр, и в 2D и 3D расчётах, а вернее при написании софта, который проводит вычисления для их решения одними из наиболее часто выполняемых операций являются векторно-матричные преобразования типа $aX+bY$, где $a,b$ — скалярные значения, $X, Y\in R^n$ — вектора или матрицы размерности $R^{n\times m}$.


Собственно вот такие:


image
(источник).


Так, чтобы не углубляться в теорию оптимизации за примерами достаточно вспомнить формулу численного интегрирования Рунге-Кутты четвёртого порядка:


$Y_{n+1}=Y_n+\frac{h}{6}(k_1 + 2 k_2 + 2 k_3+k_4),$


где $Y_i$ — очередное значение интегрируемой функции $f(t,Y)$ $h$ — шаг метода, а $k_i$, $i=1..4$ — значения интегрируемой функции в некоторых промежуточных точках — в общем случае векторах.


Как можно заметить основную массу математических операций как для векторов, так и для матриц составляют:


  • сложение и вычитание — более быстрые;
  • умножение и деление — более медленные.

О сложности вычислений хорошо написано в соответствующем курсе МФТИ.


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


Соответственно есть смысл заняться снижением количества операций привносящих наибольшую сложность — умножения (математика) и операции управления памятью (алгоритмика).


Снижение вычислительной сложности


Произвольный вектор $x\in R^n$ представим в виде $x = aX$, где
$a\in R^1$, $X\in R^n$, обозначим как $[a,X]$. Операцией приведения такого вектора назовем вычисление $[a,X]=x\in R^n$, — обратную операцию, вычислительная сложность которой составляет $n$ операций действительного умножения.


Соответственно под экономией будем понимать экономию операций умножения.


На чём можно сэкономить при операциях с векторами :


  1. Умножение на константу $b[a,X]= [ab,X]$ — экономится $n-1$ операция умножения.
  2. Сложение $k$ сонаправленных векторов $\sum\limits_{i=1}^{k} [a_i,X]= \left[\sum\limits_{i=1}^{k}a_i,X\right]$
    • откладывается вычисление $n$ умножений;
    • экономится $(k-2)n$ умножений.
  3. Сложение не сонаправленных векторов $[a,X]+[b,Y]= [a,X+b/a\ Y]$ — откладывается вычисление $n-1$ операция умножения, при $a \ne 0$.
  4. Умножение матрицы на вектор: $A[a,x]=[a,Ax]$, где $A\in R^{m\times n}$ — что позволяет отложить вычисление $m$ операций умножения.
  5. Скалярное произведение векторов $([a,X],[b,Y])= ab(X,Y)$ — экономится $n-1$ операция действительного умножения.

Аналогично, если представить матрицу в виде $Y = a\bar Y$, где
$a\in R^1$, а $\bar Y, Y\in R^{n\times m}$, обозначим как $[a,\bar Y]$. Соответственно операцией приведения — вычисление $[a,\bar Y]=Y\in R^{n\times m}$, — обратную операцию, вычислительная сложность которой составляет $n\cdot m$ операций действительного умножения.


Аналогичны и возможности экономии (1-4) с той той лишь разницей, что операция скалярного произведения (5) отсутствует, а для остальных операций возможности сэкономить умножения становится больше в $m$ раз, что довольно приятно, особенно при необходимости хранения в матричном представлении массивов векторов для проведения над ними различных массовых операций.


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


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


Снижение алгоритмической сложности


Как видно из предшествующих рассуждений, — реальная экономия возникает для:


  • умножения вектора/матрицы на число и скалярных произведений векторов;
  • операций над векторами и матрицами, отличающимися только коэффициентами.

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


class Vector {
public:
   sVec *v; //массив для данных
   double upd; //коэфф. a
   bool updated; //признак  того, что a==1
//....
};

Возможно не самым очевидное здесь — назначение переменной updated, ввиду того, что если upd = 1, то вроде как дополнительные проверки избыточны. Но если вспомнить, что последовательное деление и умножение следующего вида:


a_ /= b_;
a_ *= b_;

не обязано сохранять значение a_ из-за округлений, то видно, что данная страховка не лишняя.


Далее, если обратить внимание на, то, что если у векторов $[a,X]$ и $[b,X]$ сам $X$ одинаков, то нет смысла хранить его дважды, а достаточно создать единственный объект класса базовый вектор $X$ и ссылаться на него в расчётах, т.е.


class sVec {
public:
    unsigned long size;
    double *v;//собственно массив данных
    long linkCount;//счётчик ссылок
//..
};

image


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


Особенности реализации


Для примера (который можно взять на github) рассмотрим достаточно простые вычисления, связанные с векторно-матричной арифметикой в разрезе экономии операций умножения. Итак: сначала создадим необходимые объекты:


    Vector X(2), Y(2), Z(2); //три вектора из R^2
    double a = 2.0; // коэфф. для проверки умножения на скаляр

и инициализируем их как единичный вектор для X=[1,1] и нулевой для Y=[0,0]


    X.one(); //единичный вектор [1,1]
    Y.zero();// создаём нулевой вектор [0,0]
    cout<<"Вектор X: "<<X<<" Счётчик ссылок X: "<<X.v->linkCount<<" Адрес массива: "<<X.v<<endl;
    cout<<"Вектор Y: "<<Y<<" Счётчик ссылок Y: "<<Y.v->linkCount<<" Адрес массива: "<<Y.v<<endl;

Получаем, что внутренние вектора указывают на разыве массивы и оба счётчика ссылок равны единице.


Вектор X: [1,1]; Счётчик ссылок X: 1 Адрес массива: 448c910
Вектор Y: [0,0]; Счётчик ссылок Y: 1 Адрес массива: 448c940

Далее приравниваем вектора:


    Y=X;// теперь приравниваем X и Y
    cout<<"Вектор X: "<<X<<" Счётчик ссылок X: "<<X.v->linkCount<<" Адрес массива: "<<X.v<<endl;
    cout<<"Вектор Y: "<<Y<<" Счётчик ссылок Y: "<<Y.v->linkCount<<" Адрес массива: "<<Y.v<<endl;

и видим, что теперь оба вектора указывают на один и тот же адрес, а счётчик ссылок увеличился на 1:


Вектор X: [1,1];  Счётчик ссылок X: 2 Адрес массива: 448c910
Вектор Y: [1,1];  Счётчик ссылок Y: 2 Адрес массива: 448c910

После чего домножаем вектор X на 2


    X*=a;  //домножаем вектор на коэффициент 2.0
    cout<<"Вектор X: "<<X<<" Счётчик ссылок X: "<<X.v->linkCount<<" Адрес массива: "<<X.v<<" Коэфф. a вектора X: "<<X.upd<<endl;
    cout<<"Вектор Y: "<<Y<<" Счётчик ссылок Y: "<<Y.v->linkCount<<" Адрес массива: "<<Y.v<<" Коэфф. a вектора Y: "<<Y.upd<<endl;

и получаем, что несмотря на то, что оба вектора как и раньше указывают на один и тот же адрес, и счётчик ссылок как и прежде 2, у .X сохранённый коэффициент 2, а у Y соответственно равен 1:


Вектор X: [2,2]; Счётчик ссылок X: 2 Адрес массива: 448c910 Коэфф. a вектора X: 2
Вектор Y: [1,1]; Счётчик ссылок Y: 1 Адрес массива: 448c910 Коэфф. a вектора Y: 1

И в конце, сложив X+Y


    Z= X+Y;
    cout<<"Вектор Z: "<<Z<<" Счётчик ссылок Z: "<<Z.v->linkCount<<" Адрес массива: "<<Z.v<<" Коэфф. a вектора Z: "<<Z.upd<<endl;

получаем Z=[3,3]:


Вектор Z: [3,3]; Счётчик ссылок Z: 1 Адрес массива: 448c988 Коэфф. a вектора Z: 1

Для матриц все вычисления аналогичны:


    Matrix X_(2), Y_(2), Z_(2); //создаём 3 единичных диагональных матрицы 2x2

    cout<<"Матрица X: "<<endl<<X_<<" Счётчик ссылок X: "<<X_.v->linkCount<<" Адрес массива: "<<X_.v<<endl;
    cout<<"Матрица Y: "<<endl<<Y_<<" Счётчик ссылок Y: "<<Y_.v->linkCount<<" Адрес массива: "<<Y_.v<<endl;

    Y_=X_;// теперь приравниваем X_ и Y_
    cout<<"Матрица X: "<<endl<<X_<<" Счётчик ссылок X: "<<X_.v->linkCount<<" Адрес массива: "<<X_.v<<endl;
    cout<<"Матрица Y: "<<endl<<Y_<<" Счётчик ссылок Y: "<<Y_.v->linkCount<<" Адрес массива: "<<Y_.v<<endl;

    X_*=a;  //домножаем матрицу на коэффициент 2.0
    cout<<"Матрица X: "<<endl<<X_<<" Счётчик ссылок X: "<<X_.v->linkCount<<" Адрес массива: "<<X_.v<<" Коэфф. a матрицы X: "<<X_.upd<<endl;

    Z_= X_+Y_;
    cout<<"Матрица Z: "<<endl<<Z_<<" Счётчик ссылок Z: "<<Z_.v->linkCount<<" Адрес массива: "<<Z_.v<<" Коэфф. a матрицы Z: "<<Z_.upd<<endl;

Вывод также аналогичен:


Матрица X:
[1,0]
[0,1];
 Счётчик ссылок X: 1 Адрес массива: 44854a0
Матрица Y:
[1,0]
[0,1];
 Счётчик ссылок Y: 1 Адрес массива: 44854c0
Матрица X:
[1,0]
[0,1];
 Счётчик ссылок X: 2 Адрес массива: 44854a0
Матрица Y:
[1,0]
[0,1];
 Счётчик ссылок Y: 2 Адрес массива: 44854a0
Матрица X:
[2,0]
[0,2];
 Счётчик ссылок X: 2 Адрес массива: 44854a0 Коэфф. a матрицы X: 2
Матрица Z:
[3,0]
[0,3];
 Счётчик ссылок Z: 1 Адрес массива: 4485500 Коэфф. a матрицы Z: 1

Библиотечка доступна на GitHub.


Пользуйтесь, критикуйте.

Tags:
Hubs:
-1
Comments 13
Comments Comments 13

Articles