Подход к решению систем линейных алгебраических уравнений с интервальной неопределенностью в исходных данных

Система, линейных алгебраических уравнений - это фундаментальный объект, который встречается при решении многих задач. Часто оказывается, что коэффициенты рассматриваемой системы не могут быть заданы точно, но известны интервалы, которым они принадлежат. В условиях такой интервальной неопределенности коэффициентов необходимо уточнение определения решения. В дальнейшем изложении будем использовать стандартную нотацию, принятую в интервальном анализе [1].

В работах [1-5] систематизированы подходы к учету интервальной неопределенности и дана их классификация. В соответствии с данной классификацией, AS -решением системы линейных алгебраических уравнений A x = b, в которой элементы матриц A и b представляют интервалы aij = [ a ^ , a ^ ], bj = j , ,b j^, i,j = 1 , 2 ,... ,n, называют элементы допускового множества.

n

0 toi (A , b) = x : ( Vi, j = 1 , 2 ,...,n ) ( Va j G aij)        a ij X j G bi

\ j =1

SS -решением рассматриваемой системы уравнений называют точки объединенного множества.

n

⁰ uni ^(A , ^b) — * ^x :  ^{( Vi}, j — ¹ , ² , . . . , n ) ⁽ ^^a ij G aij) [       ^a ij x j G  ^bi

\ j =1

В работах [6, 7] доказано, что поиск SS -решения интервальной системы линейных уравнений является NP-трудной задачей. С другой стороны, в соответствии с теоремой Рона [8] любая точка допустимого множества AS -решений допускает представление в виде x = x + — x . где x + , x^- являются решением системы неравенств

n

Z [aij x+ j=1

^— a ij x j   > b i ,

n

^^ a iii j x+ — a ij x— < b i , i = 1 , 2 ,... ,n, x + , x- >  0 . j =1

Следовательно, задача поиска AS -решений имеет полиномиальную сложность.

Методы оценки AS -решений для случаев 0 toi ( A , b ) = 0 рассмотрены в работах [15] и др. Основным методом исследования допускового множества решений, развиваемым в Новосибирске [3-5], является « метод распознающего функционала » . В нем для принятия решения о разрешимости или неразрешимости задачи (т.е. о пустоте/непустоте множества решений) необходимо поработать с некоторым специальным (негладким и вогнутым) функционалом, который назван « распознающим » . При этом максимизация распознающего функционала, которую практически можно выполнять, например, с помощью различных методов негладкой оптимизации, разработанных в Институте кибернетики НАН Украины [9], дает достаточно содержательную информацию для возможной коррекции задачи. Разработанные С.П. Шарым и П.И. Стецюком программы для исследования разрешимости интервальной линейной задачи о допусках (пустоты-непустоты допускового множества решений), имеются в свободном доступе на сайте [10]. Программы реализованы в INTLAB’e - интервальном расширении MATLAB’a, а также в Int4Sci - интервальном расширении Scilab’a.

Во многих практических задачах система неравенств (1) оказывается плохо обусловленной или вообще несовместной. В этом случае по аналогии с работами [11, 12] разумным представляется введение понятия « псевдорешения ».

Целью данной работы является изложение анонсированного в работе [13] понятия « псевдорешение » для систем уравнений с интервальной неопределенностью и способы построения инструментальных программных средств их поиска.

1.    Псевдорешение системы интервальных уравнений

Пусть дана система линейных алгебраических уравнений A x = b, в которой элементы матриц A и b представляют интервалы aij = [ aij, aij ]. bj = bb, ,bj]. i,j = 1 , 2 ,...,n.

Для заданной системы уравнений построим параметризованное семейство систем уравнений A x = b( z ) c модифицирован!юй правой частью b( z ) = [ b — z |b|, b + z | b |] ,z >  0.

Пусть z* = inf {z : 0 toi (A , b( z )) = 0} . Псевдорешеиием исходной системы A x = b будем называть внутренние точки допустимого множества 0 toi (A , b( z* )).

Корректность введенного определения подтверждает

Теорема 1. Для любой системы интервальных уравнений A x = b при всех z >  1 мномсе-ство 0 toi ( A , b ( z )) = 0.

Доказательство. В соответствии с теоремой Рона условие 0 toi (A , b( z )) = 0 эквивалентно совместности системы линейных неравенств

n

^ | a ij x + — a ij x— > b i — z \b i \, i = 1 , 2 ,...,n., j =1

n

^ [ « у x j — » ij x-  <  bt + z\b i \, i = 1 , 2 ,...,n,                   (2)

j =1

x ⁺ , x^-  >   0 .                                                 (3)

Полагая в (1) - (3) x ⁺ = x = 0. получим

0 > b i — zlb i |, 0 < b t + zb i |, i = 1 , 2 , ...,n.

Таким образом, для всех z >  1 имеет место включение 0 € 0 toi ( A , b ( z )). Теорема доказана.

Т

• 2.    Способ поиска псевдорешения

Способ нахождения псевдорешения системы уравнений A x = b дает

Теорема 2. Существует решение x ⁺ , вания

z

n

52⁽ a tj ^x+ ^— a tj ^x- )

j =1

n

52 ⁽ a ij ^x+ ^— a ij x j )

j =1

^{x +} ,^x-,z

x  € R ⁿ, z* € R задачи линейного программиро-

^ min ,(4)

x ⁺ , x ^- , z

• >  bi — ^ЬД i = 1, 2 ,...,п,(5)

• <   bi + z|bi |, i = 1, 2 ,...,n,(C)

• >  0, j = 1, 2 ,...,n,(7)

3.    Техника реализации

при этом x* = x ⁺ — x является псевдорешением системы A x = b .

Доказательство. Сначала докажем существование оптимального решения x⁺ , x^- € R ⁿ, z* € R задачи линейного пограммирования (4)-(7). Из теоремы 1 и теоремы Рона следует, что множество допустимых решений рассматриваемой задачи не пусто. Задача, двойственная рассматриваемой, имеет вид

n 52 bty 11" i=1 n 52 biy 2 i i=1 → max, y1i^y2i (8) n 52 j 11 — n 52 aj^y 2 i ≤ 0, j = 1, 2 ,...,n, (9) i=1 n — ^ajiy 1 i + i=1 n 52 ajiy 2 i ≤ 0, j = 1, 2 ,...,n, (Ю) i=1 n ^ |bi |y 1i + i=1 n 52 |bi|y2i ≤ 1, i = 1, 2,... ,n, (П) i=1 i=1 y 1 i,y 2 i ≥ 0, i = 1, 2 ,...,n. (12) Легко заметить, что решение y 1 i = ; y 2 i = 0, i = 1, 2,... ,n является допустимым решением задачи (8)—(12). Таким образом, показано существование допустимых решений как у прямой, так и двойственной задач линейного программирования (задачи (4)-(7) и (8)—(12)). Из теоремы двойственности в линейном программировании следует существование у этих задач оптимальных решений.

Пусть x + , x- , z* оптимальное решение задачи (4)-(7). Из теоремы Рона следует, что x* = x + — x- является допустимым решением системы A x = b ( z* ). Из оптимальности z* следует, что x* является псевдорешением интервальной системы линейных уравнений A x = b Теорема доказана. □

Таким образом, введенное понятие « псевдорешение » интервальной системы линейных уравнений является вполне конструктивным и позволяет давать результаты при решении интервальных систем в том случае, когда пусковое множество 0 toi ( A , b ) допустимых решений пусто.

Однако, следует обратить внимание на высокую степень вырожденност задач (4)-(7) и (8)—(12)), что будет приводить при использовании приближенных вычислений к зацикливанию симплекс-метода. Избежать зацикливания можно за счет использования вычислений без округления [14, 15]. В этом случае на каждой итерации симплекс-метода количество требуемых бит памяти не превосходит величины 4 lm ⁴ + O ( lm 3), г де m - минимальная из размерностей задачи, l - число бит, достаточных для представления одного элемента матрицы исходных данных, при этом эффективность распараллеливания (т.е. отношение ускорения к числу процессоров) составляет в асимптотике величину, близкую к 100% [16].

Итак, для успешной реализации изложенного способа поиска псевдорешения на компьютере понадобится несколько составляющих. С одной стороны необходимо обеспечить достаточную точность вычислений для преодоления зацикливания симплекс-метода. Подобная техника описана в [17]. Далее, для обеспечения эффективности всего алгоритма в целом, следует обеспечить достаточное быстродействие точных вычислений, а для использования параллельной архитектуры современных процессоров, произвести эффективную декомпозицию задачи линейного программирования на потоки.

• 3.1.    Обеспечение необходимой точности

  В рамках предыдущих исследований были созданы классы overlong и rational, реализованные в объектно-ориентированной парадигме на языке С — как библиотека классов Exact Computation [14]. Данные классы позволяют производить безошибочные дробнорациональные вычисления и имеют следующие характеристики. Объектами класса rational являются обыкновенные дроби p/q, где p, q - объекты класса overlong. Класс overlong предназначен для расширения логических возможностей целочисленных вычислений на компьютере. Объем памяти, занимаемый такими объектами, определяется значениями представляемых чисел, их диапазон ограничен только объемом адресуемой памяти. Диапазон чисел, представляемых объектами класса overlong, расширен до (2 - ³² * ⁶⁵⁵³⁵, 2³² * 6⁵⁵³⁵), _а минимальный шаг дискретизации чисел, представляемых объектами класса rational, может достигать (2 - ^2097150), Для объектов классов overlong и rational определены все операторы, операции и бинарные отношения, используемые для стандартных числовых типов данных. Таким образом, классы overlong и rational дают потенциальную возможность использовать в программах пользователя безошибочное выполнение основных арифметических операций над полем рациональных чисел.

На сегодняшний день возможность использования безошибочных вычислений представляет известная библиотека GMP(The GNU Multiple Precision Arithmetic Library) [15]. Библиотека распространяется под лицензией GNU LGPL, актуальная версия библиотеки GMP 5.1.1 доступна для загрузки с официального сайта проекта. Программный код оптимизирован под большинство существующих процессорных архитектур, однако она не предо- ставляет своим объектам возможность их использования в распределенных вычислениях. Для более полноценного использования современных процессорных архитектур классы overlong и rational хранят и оперируют числами по основанию 232, код операций оптимизирован для системы счисления по основанию степень двойки. Оптимизации применяются также при работе с памятью, поскольку в СН—h нет автоматического сборщика мусора, то лишние перевыделения памяти приводят ее фрагментации и снижению быстродействия приложения в целом, краткое описание современных реализаций классов дано в [18]. Для облегчения сопровождения и модификации классов операции с памятью инкапсулированы в отдельный класс MemHandle, а выполнение базовых арифметических операций с данными полностью производится в рамках класса Arif Realization (см. фрагмент листинга 1). Тем

class overlongNM { private :

static Ar if Realization realization;

private :

MemHandle mhandle ;

.   .   .

public: inline int32 size() const {return leng;}//leng of number public: inline int32 sign() const {return sgn;} //sign of number e ■   ■

//addition templatefriend const overlongNM operator -

{overlongNM rez(num); return (rez+=v);} friend const overlongNM operatorЬ

(const overlongNM&, const overlongNM&);

■   ■   ■

}

Листинг 1. Фрагмент класса overlong самым объект класса overlong содержит в себе объект типа MemHandle, и все действия с памятью происходят через интерфейс MemHandle. Все арифметические операции с данными осуществляются вызовом соответствующих методов класса Arif Realization. Пример реализации метода add класса overlong представлен в листинге 2. Эта техника позволяет пол-

void overlongNM :: add ( const overlongNM& alpha , const overlongNM& beta){ d_t carry;

const overlongNM& a=( alpha . s ize ()> = bet a . s ize () ) ? alpha : beta ;

const overlongNM& b=( alpha . s ize ()> = bet a . s ize () ) ? beta : alpha ;

//вызов базовой арифметической операции

ArifRealization : : add ( a.. mhandle . getptr ( ) ,LA, b . mhandle . getptr () ,LB, mhandle . providetmpptr (LA, 1) , newleng , carry ) ;

mhandle . settmpasptr ( ) ;

if ( carry ) mhandle . safeset value (LA, carry ) ;

leng=newleng ;

sgn=sg;

}

Листинг 2. Организация операции сложения ностью абстрагироваться от реального места хранения данных (за это полностью отвечает реализация класса MemHandle) и способа реализации арифметических операций (конкретная методика выполнения операций над разрядами инкапсулирована в Arif Realization).

• 3.2.    Обеспечение необходимой производительности

• 3.3.    Мелкозернистый параллелизм

• 3.3.1.    Сложение

• 3.3.2.    Умножение

Вышеописанная техника разбиения класса overlong, на три части «интерфейс-память-арифметика», позволяет гибко использовать возможности вычислительной системы. Поскольку производительность алгоритма симплекс метода обеспечивается, в том числе, за счет эффективности реализаций классов overlong и rational, при их написании учтена возможность использования современных гетерогенных вычислительных сред с GPU ускорителями Nvidia (возможность задействовать GPU от компании AMD и встроенные графических ускорителей, например, Intel HD Graphics дает язык OpenCL, однако эксперименты показали большую трудоемкость процесса кодирования при малом приросте эффективности). Алгоритмы параллельного выполнения базовых арифметических операций, а также некоторые аспекты их реализации в гетерогенной среде описаны в [18]. Хранение операндов организуется с учетом устройства, на котором производятся вычисления. Так, при наличии в системе GPU Nvidia, все данные (разряды) чисел можно хранить непосредственно на GPU, там же выполнять арифметические операции над числами, это снижает до минимума количество пересылок данных по PCI шине, остается лишь скопировать в основную память результаты вычислений. Для систем без GPU, вычисления проводятся на процессоре, данные хранятся непосредственно в оперативной памяти системы.

Параллельные вычисления на GPU требуют переработки алгоритмов базовых арифметических операций с учетом специфической архитектуры устройства. Ниже приводятся листинги некоторых операций для GPU Nvidia на расширениии языка С от Nvidia Corporation (CUDA С).

Операция сложения длинных чисел на GPU осуществляется в несколько этапов: параллельное сложение разрядов, синхронизация, параллельное распространение переносов из разрядов. Особенность архитектуры GPU, а именно, выполнение нитей блоками и отсутствие синхронизации между блоками требует сохранения « пограничных » переносов во временом массиве (за это отвечает строка номер 16 в листинге 3 и последующего распространения (DNumAdd_part2) (листинг 4). Установка параметров и запуск кода на графическом ускорителе представлены на листинге 5. Этот код выполняется на стороне CPU, или, так называемой, стороне Host.

Операция умножения одна из наиболее затратных по времени. Для выполнения операции параллельного умножения используется быстрая разделяемая между нитямя блока __ shared __память. В реализации существенно используется особенности архитектуры GPU Nvidia, а именно, полностью синхронное выполнение инструкций в рамках одного warp. Это избавляет от необходимости выполнять синхронизацию между нитями. Исходный код для выполнения на GPU (kernel) представлен на листинге 6. Окончательное формирование результата происходит последовательным проходом по массиву rez[] и преобразованием 64-битного числа rez[j] в непосредственно разряд ответа и разряд переноса. Эти действия

__global__void DNumAdd_partl

(d_t * A. int32 LA.d_t * B.int32 LB, d_t * C.d_t * bGCarry . int32 * f ) { int 3 2 gId=blockDim . x * blockldx.x + threadldx.x; int 64 tmp = 0;

if (gid >= LA) return; //bound check if(g!d>=LB) C[ g!d]=A[ gid ] ;

else {tmp=(int64 )A[ gid] ■ (int64 )B [ gid ] ; C[ gId]=tmp&MAX_DIGIT;} _syncthreads ();// carry propagation in the block int 3 2 Hd=threadldx . x ■ 1, i=gld ■ 1, gS=blockDim . x ;

for (tmp»=BIT_IN_DIGIT; tmp 11 d -4 , i -4 ){ if (1 Id==gS ) { bGCarry [ blockldx . x]=tmp ; *f = l; return;} tmp+=( int 64 )C [ i ]; C[ i]=tmp&MAX_DIGIT; tmp»=BIT_IN_DIGIT;

} i f ( i==LA & tmp ) {bGCarry [ Lcarry ] = 1; }

}

Листинг 3. Поразрядное сложение

__global__void DNumAdd_part2

(d_t * C, d_t * bGCarry, int 3 2 gS , uint32 Lcarry, uint32 LA) {

//carry propagation between blocks int gid = blockDim . x *blockldx.x + threadldx.x, i=(gId + l)*gS;

if (gid >= Lcarry) return;

uint64 tmp=(uint64 ) bGCarry [ gid ] ;

for (; tmp & i

C [ i ] = tmp&MAX_DIGIT;

tmp»=BIT_IN_ DIGIT;

} if ( i==LA & tmp ) {bGCarry [ Lcarry ] = 1; }

}

Листинг 4. Параллельное распространение переносов также выполняются на GPU, но в однопоточном режиме.

• 3.3.3.    Деление

• 3.4.    Крупнозернистый параллелизм

4.    Вычислительный эксперимент

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

Применяемая техника параллельной реализации симплекс метода детально описана в [19]. В данном случае для разбиения задачи на потоки использовался механизм многопоточного програмирования std::thread, предоставляемый стандартом С 44 0x11.

Вычислительный эксперимент проводился на компьютере с процессором Intel Core i7-950 3.06 ГГц, 6 Гб ОЗУ, GPU Nvidia 460(1гб GDDR5), под управлением ОС Win 7 х64, в качестве компилятора был выбран 64-разрядный Visual С 1—4 2011. В качестве модельной

void ArifRealization :: add ( const d_t *A, int32 LA, const d_t *B,int32 LB,d_t *C, int32 d_t *bGCaiTy_d=ATLL. * ansCarry_li=new d_t [1] :

i ii 13 2 * cF _ d=XULL. * cF _ li=new hi 13 2 [ 1 ] :

cudaMalloc (( void **)УЬСCarry_d . s iz c of (d_t) * ( bPerGrid 1));

cudaMemset (( void *)bGCarry_d . 0. s iz c of (d_t) * ( bPerGrid 1));

cudaMalloc (( void **)&cF_d , sizeof(int32));

cudaMemset (( void *)cF_d, 0, s iz e of (int 3 2 ));

DNumAdd_partl «< bPerGrid , tPerBlock »>

( const _cast *>(A) , LA, const_cast*>(B) , LB, C, bGCarry_d ,CF_d) ;

cudaMemcpy (cF_h, cF_d , sizeof (int32 ) , cudaMemcpyDeviceToHost ) ;

if( *cF_h){'                                         '

int gS=tPerBlock , LCarry=bPerGrid ;

bPerGrid = (bPerGrid ‘tPerBlock — 1) tPerBlock:

DNumAdd_part2 «< bPerGrid , tPerBlock:-»

(d_buffC, bGC_d, gS , LCarry, LA);

} cudaMemcpy ( ansCarry_h ,&bGCarry_d [ bPerGrid ] , sizeof (d_t), cudaMemcpyDeviceToHost ) ;

ND= (carry = *ansCarry_h) ? LA 1: LA;

delete [] ansCarry_h; cudaFree ( bGCarry_d ); cudaFree(cF_d );

}

Листинг 5. Вызов функции для исполнения на GPU

__global__void DNumMult(d_t *A,int32 LA,d_t *B,int32 LB,d_t *rez) { int32 Hd=threadldx . x , gId=blockDim . x *blockldx.x + lid;

if(gId>=LB) return;

int 32 cBS=(LA blockDim . x — 1)/blockDim . x ;

_shared__uint64 sha [] , shrez [ ] ;

for (int i = lld *cBS:i (lid l)*cBS MC 1 LA: 1   )( s h a [ i ] = A [ i ] ; s h r e z [ i ] = 0;

} shrez [LA lid [ = 0;

uint64 d igit =( uint64)B[gid] , t=0UL;

fоr(i n t i = 0; i

shrez [ i + Hd] k=t£MAX_DIGIT;

t»=BIT_IN_DIGIT;

} shrez [LA lid] -=t cBS = (LA blockDim . x blockDim . x — l)/blockDim . x ;

for (int i = lld *cBS;i<(Hd l)*cBS && i

Atomic Add ( rez [ i+g Id ] , shrez [ i ]) ;

}

}

Листинг 6. Умножение на GPU задачи использована система с матрицами:

A =

Г i * (1 — 5) i * (1 + 5) 1

I i + j —¹, i + j —¹ Lxn

; b = [1, 1 /2,..., 1 /(n — 1), 1 /n]T.

Зависимость минимального расширения правой части (параметр z*), соответствующего псевдорешению, при фиксированном значении n = 20 приведена в таблице 1.

Таблица 1

Минимальное расширение правой части системы

δ	10-1	10-2	10-3	10-4	10-5	10-6
z∗	0,81	0,389	ОД	0,025	0,0062	0,0017

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

Таблица 2

Время работы

Размерность матрицы (и)	10	20	50	100
Время работы	0,46 с	7,73 с	7,39 м	15,1 ч

5.    Заключение

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

Исследование выполнено при поддермске Министерства образования и науки Российской Федерации, соглашение Ц.631.21.0395.