Лек­ция № 4: При­ве­де­ние мат­ри­цы к тре­уголь­но­му ви­ду

При­ве­де­ние мат­ри­цы к тре­уголь­но­му ви­ду тре­бу­ет­ся, пре­жде все­го, для ре­ше­ния си­стем ли­ней­ных урав­не­ний (СЛАУ). Од­на­ко, для про­сто­ты весь ма­те­ри­ал лек­ции бу­дет рас­смот­рен на при­ме­ре вы­чис­ле­ния опре­де­ли­те­ля мат­ри­цы. Ведь са­мый быст­рый спо­соб вы­чис­лить опре­де­ли­тель — при­ве­сти мат­ри­цу к тре­уголь­но­му ви­ду, а за­тем пе­ре­мно­жить её диа­го­наль­ные эле­мен­ты. Ес­ли же вы хо­ти­те при­ме­нить ма­те­ри­ал этой лек­ции к ре­ше­нию СЛАУ, то един­ствен­ное что вам нуж­но сде­лать — до­ба­вить к опе­ра­ци­ям со стро­ка­ми мат­ри­цы та­кие же опе­ра­ции с эле­мен­та­ми пра­вой ча­сти СЛАУ (ну и не за­быть най­ти зна­че­ния не­из­ве­ст­ных ве­ли­чин с по­мо­щью об­рат­но­го хо­да ал­го­рит­ма Гаус­са).

Вез­де в этой лек­ции стро́ки и столб­цы мат­риц ну­ме­ру­ют­ся с ну­ля, что­бы бо­лее со­от­вет­ство­вать язы­ку про­грам­ми­ро­ва­ния C++. Эле­мен­ты мат­ри­цы обо­зна­ча­ют­ся, как a_{ij}, где i=0,...,N-1 — но­мер строки́, j=0,...,N-1 — но­мер столб­ца. Для со­кра­ще­ния за­пи­си стро́ки мат­ри­цы обо­зна­ча­ют­ся, как \vec a_i: с ни­ми мож­но про­из­во­дить опе­ра­ции так же, как с век­то­ра­ми.

Ал­го­ритм Гаус­са

Мат­ри­цу мож­но при­ве­сти к тре­уголь­но­му ви­ду при по­мо­щи ал­го­рит­ма Гаус­са. Идея ал­го­рит­ма за­клю­ча­ет­ся в том, что эле­мен­ты a_{ij}, i>j об­ну­ля­ют­ся пу­тём вы­чи­та­ния строки́ но­мер j, умно­жен­ной на a_{ij}/a_{jj}, из строки́ но­мер i. Фор­маль­но ал­го­ритм вы­гля­дит так:

\begin{array}{l}
\text{For }j=0,...,N-2\\
\left|\quad\begin{array}{l}
	\text{For }i=j+1,...,N-1\\
	\left|\quad
		\vec a_i \leftarrow \vec a_i-\cfrac{a_{ij}}{a_{jj}}\cdot\vec a_j
	\right.\\
\end{array}\right.\\
\end{array}
(1)

Ал­го­ритм стро­ит верх­не­тре­уголь­ную мат­ри­цу, из­ме­няя ис­ход­ную мат­ри­цу. Ес­ли ис­ход­ная мат­ри­ца долж­на со­хра­нить­ся, то пе­ред на­ча­лом ра­бо­ты ал­го­рит­ма нуж­но со­здать её ко­пию.

Не­смот­ря на то, что в ал­го­рит­ме (1) по­ка­за­ны два вло­жен­ных цик­ла, на са­мом де­ле в нём три цик­ла, так как вы­чи­та­ние строк — опе­ра­ция над не­сколь­ки­ми эле­мен­та­ми. Со­от­вет­ствен­но, вре­мя ра­бо­ты ал­го­рит­ма со­став­ля­ет O\left(N^3\right).

Ал­го­ритм Гаус­са не бу­дет ра­бо­тать, ес­ли в ка­кой-то мо­мент ока­жет­ся a_{jj}=0 в фор­му­ле (1). Не­об­хо­ди­мое и до­ста­точ­ное усло­вие то­го, что в ал­го­рит­ме не воз­ник­нет де­ле­ния на ноль — все глав­ные ми­но­ры мат­ри­цы долж­ны быть от­лич­ны от ну­ля.

По­сле при­ве­де­ния мат­ри­цы к верх­не­тре­уголь­но­му ви­ду её опре­де­ли­тель мож­но вы­чис­лить пе­ре­мно­же­ни­ем диа­го­наль­ных эле­мен­тов:

\begin{array}{l}
d\leftarrow a_{00}\\
\text{For }i=1,...,N-1\\
\left|\quad
	d\leftarrow d\cdot a_{ii}
\right.\\
\end{array}
(2)

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

При­мер ра­бо­ты ал­го­рит­ма Гаус­са

Вы­чис­лим опре­де­ли­тель сле­дую­щей мат­ри­цы по при­ве­дён­но­му вы­ше ал­го­рит­му:

\begin{pmatrix}3&2&3&4\\
4&4&3&2\\
1&4&4&3\\
2&3&1&1\end{pmatrix}
(3)

Пред­по­ло­жим, что ком­пью­тер при вы­чис­ле­ни­ях хра­нит не бо­лее 5-ти зна­ча­щих цифр: так бо­лее на­гляд­но вид­ны на­кап­ли­ваю­щие­ся ошиб­ки округ­ле­ния. Умно­жим ну­ле­вую стро­ку на 4/3 и вы­чтем её из пер­вой строки́:

\begin{pmatrix}3&2&3&4\\
0&1.3333&-1&-3.3333\\
1&4&4&3\\
2&3&1&1\end{pmatrix}
(4)

Ана­ло­гич­но вы­чтем ну­ле­вую стро­ку из вто­рой и тре­тьей строк (умно­жив её пе­ред этим на 1/3 и 2/3 со­от­вет­ствен­но):

\begin{pmatrix}3&2&3&4\\
0&1.3333&-1&-3.3333\\
0&3.3333&3&1.6667\\
0&1.6667&-1&-1.6667\end{pmatrix}
(5)

Вы­чтем первую стро­ку из вто­рой и тре­тьей (стро́ки у нас ну­ме­ру­ют­ся с ну­ля!), умно­жив её пе­ред этим на 3.3333/1.3333 и 1.6667/1.3333 со­от­вет­ствен­но:

\begin{pmatrix}3&2&3&4\\
0&1.3333&-1&-3.3333\\
0&0&5.5&10.0001\\
0&0&0.2501&2.5001\end{pmatrix}
(6)

На­ко­нец, вы­чтем из тре­тьей строки́ вто­рую, умно­жен­ную на 0.2501/5.5:

\begin{pmatrix}3&2&3&4\\
0&1.3333&-1&-3.3333\\
0&0&5.5&10.0001\\

0&0&0&2.0454\end{pmatrix}
(7)

Пе­ре­мно­жая диа­го­наль­ные эле­мен­ты мат­ри­цы, по­лу­чим зна­че­ние опре­де­ли­те­ля: 44.9977. Не­труд­но до­га­дать­ся, что ес­ли бы мы вы­чис­ля­ли опре­де­ли­тель точ­но (на­при­мер, с ис­поль­зо­ва­ни­ем дро­бей), то он был бы ра­вен 45.

Пе­ре­ста­нов­ка строк

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

За­ме­чу, что ну­ле­вой зна­ме­натель в вы­ра­же­нии (1) ещё не озна­ча­ет ра­вен­ства ну­лю опре­де­ли­те­ля всей мат­ри­цы; так­же этот факт не озна­ча­ет не­воз­мож­но­сти ре­ше­ния СЛАУ с дан­ной мат­ри­цей. Кро­ме то­го, воз­мож­на та­кая си­ту­а­ция, ко­гда зна­ме­натель в вы­ра­же­нии (1) из-за оши­бок округ­ле­ния не ра­вен ну­лю в то вре­мя, ко­гда он дол­жен быть ра­вен ну­лю. В этом слу­чае ал­го­ритм даст ре­зуль­тат, сколь угод­но силь­но от­ли­чаю­щий­ся от точ­но­го.

Для борь­бы с опи­сан­ны­ми вы­ше про­бле­ма­ми ис­поль­зу­ет­ся пе­ре­ста­нов­ка строк: пе­ред каж­дым вы­чи­та­ни­ем строки́ из всех по­сле­дую­щих (пе­ред вто­рым цик­лом ал­го­рит­ма (1)) на­хо­дит­ся стро­ка но­мер k (сре­ди строк i=j,...,N-1) с мак­си­маль­ным по мо­ду­лю эле­мен­том a_{jk}, и, ес­ли k\ne j, про­из­во­дит­ся пе­ре­ста­нов­ка строк k и j. По­сле пе­ре­ста­нов­ки стро­ка с мак­си­маль­ным по мо­ду­лю пер­вым не­ну­ле­вым эле­мен­том ока­зы­ва­ет­ся на ме­сте той строки́, ко­то­рая вы­чи­та­ет­ся из всех сле­дую­щих. И в зна­ме­натель вы­ра­же­ния (1) по­па­да­ет мак­си­маль­ный эле­мент сре­ди всех, ко­то­рые на­хо­дят­ся сни­зу от a_{jj}.

\begin{array}{l}
\text{For }j=0,...,N-2\\
\left|\quad\begin{array}{l}
	\text{Find line }k\ \left(k=j,...,N-1\right)\text{ having maximum }\left|a_{jk}\right|\\
	\text{If }k\ne j\text{ then swap lines }k\text { and }j\\
	\text{For }i=j+1,...,N-1\\
	\left|\quad\vec a_i \leftarrow \vec a_i-\cfrac{a_{ij}}{a_{jj}}\cdot\vec a_j\right.\\
	\text{End}
\end{array}\right.\\
\text{End}
\end{array}
(8)

При вы­чис­ле­нии опре­де­ли­те­ля по тре­уголь­ной мат­ри­це, сге­не­ри­ро­ван­ной ал­го­рит­мом (8), нуж­но пом­нить о том, что опре­де­ли­тель ме­ня­ет знак при каж­дой пе­ре­ста­нов­ке строк. По­это­му нуж­но по­счи­тать чис­ло пе­ре­ста­но­вок, и умно­жить про­из­ве­де­ние диа­го­наль­ных эле­мен­тов на -1, ес­ли чис­ло пе­ре­ста­но­вок не­чёт­но.

Ес­ли при­ме­нить ал­го­ритм (8) к мат­ри­це (3), то по­лу­чим сле­ду­ю­щую мат­ри­цу (с учё­том округ­ле­ния всех про­ме­жу­точ­ных ре­зуль­та­тов до пя­ти зна­ча­щих цифр):

\begin{pmatrix}4&4&3&2\\
0&3&3.25&2.5\\
0&0&1.8333&3.3333\\
0&0&0&2.0454\end{pmatrix}
(9)

За вре­мя ра­бо­ты ал­го­рит­ма бы­ло про­из­ве­де­но 2 пе­ре­ста­нов­ки, по­это­му опре­де­ли­тель ра­вен 44.9980, что не­мно­го точ­нее, чем ре­зуль­тат, по­лу­чен­ный ра­нее.

Мож­но сде­лать ал­го­ритм ещё точ­нее, до­ба­вив пе­ре­ста­нов­ку столб­цов. Од­на­ко, пе­ре­ста­нов­ка столб­цов ред­ко ис­поль­зу­ет­ся на прак­ти­ке, так как в слу­чае ре­ше­ния си­сте­мы урав­не­ний она со­про­вож­да­ет­ся пе­ре­ста­нов­кой не­из­ве­ст­ных ве­ли­чин.

Про­грамм­ная реа­ли­за­ция ал­го­рит­ма Гаус­са

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

Ст­рук­ту­ра дан­ных для хра­не­ния мат­ри­цы

Ри­су­нок 1. Ст­рук­ту­ра дан­ных для хра­не­ния мат­ри­цы

Эта ст­рук­ту­ра хо­ро­ша тем, что стро́ки мат­ри­цы мож­но по­ме­нять ме­ста­ми очень быст­ро, про­сто по­ме­няв ме­ста­ми зна­че­ния со­от­вет­ствую­щих ука­за­те­лей.

Пусть мы про­грам­ми­ру­ем на C++, и чис­ла мат­ри­цы име­ют не­ко­то­рый тип T, то­гда ука­за­те­ли на стро́ки бу­дут иметь тип T*, а са­ма мат­ри­ца M бу­дет иметь тип T**. Учи­ты­вая это, ал­го­ритм (8) мо­жет быть реа­ли­зо­ван сле­дую­щим об­ра­зом:

template<class T> T GaussDeterminant(T **const M, int const N)
{
    T determinant(1);
    int exchanges(0);
   
    for(int l1(0); l1<N-1; ++l1)
    { //Перебираю все строки матрицы, кроме последней
        int maxN( l1 );
        T maxValue( fabs( M[l1][l1] ) );
        for(int l2(l1+1); l2<N; ++l2)
        { //Нахожу строку с максимальным по модулю элементом
            T const value( fabs( M[l2][l1] ) );
            if( value > maxValue ) { maxN=l2; maxValue=value; }
        }
       
        if( maxN > l1 )
        { //Нужен обмен
            T *const temp(M[l1]); M[l1]=M[maxN]; M[maxN]=temp;
            ++exchanges;
        } else { //Обмен не нужен, но нужна проверка на ноль
            if( maxValue == T(0) ) return maxValue;
        }

        T const value( M[l1][l1] );
        determinant *= value;

        for(int l2(l1+1); l2<N; ++l2)
        { //Вычитаю строку из всех последующих
            T const k( M[l2][l1]/value );
            M[l2][l1] = T(0);
            for(int c(l1+1); c<N; ++c) M[l2][c] -= M[l1][c]*k;
        }
    }
    determinant *= M[N-1][N-1];
    if(exchanges%2) return -determinant; else return determinant;
}

Ли­стинг 1. Реа­ли­за­ция ал­го­рит­ма Гаус­са с пе­ре­ста­нов­кой строк

Мо­дифи­ка­ция Ба­рей­са

Ес­ли мат­ри­ца со­сто­ит из це­лых чи­сел, ал­го­ритм Ба­рей­са поз­во­ля­ет при­ве­сти её к тре­уголь­но­му ви­ду с ис­поль­зо­ва­ни­ем толь­ко це­ло­чис­лен­ной ариф­ме­ти­ки!

Мы ви­де­ли, что са­мый боль­шой вклад в на­кап­ли­ва­е­мую по­греш­ность да­ёт де­ле­ние на a_{ii} в ал­го­рит­мах (1) и (8). По­пы­та­ем­ся из­ба­вить­ся от это­го де­ле­ния. Пре­жде все­го за­ме­тим, что ес­ли сра­зу по­сле вы­чи­та­ния од­ной строки́ из дру­гой в фор­му­ле (1) умно­жить ре­зуль­тат на a_{ii}, то по­лу­чит­ся вы­ра­же­ние, не со­дер­жа­щее де­ле­ния:

\vec a_i \leftarrow a_{jj}\cdot\vec a_i-a_{ij}\cdot\vec a_j.
(10)

Уже с та­кой мо­дифи­ка­ци­ей ал­го­ритм мож­но ис­поль­зо­вать для ре­ше­ния СЛАУ, од­на­ко дан­ный ме­тод име­ет сле­дую­щие про­бле­мы:

  • Мо­дифи­ка­ция эк­ви­ва­лент­на мно­го­крат­но­му умно­же­нию строк мат­ри­цы на её эле­мен­ты. Ес­ли эле­мен­ты мат­ри­цы боль­ше еди­ни­цы, то для мат­риц боль­шо­го раз­ме­ра чи́сла к кон­цу ра­бо­ты ал­го­рит­ма вы­рас­тут очень зна­чи­тель­но. Воз­мож­но пе­ре­пол­не­ние ис­поль­зу­е­мо­го ти­па дан­ных. Ана­ло­гич­но, ес­ли эле­мен­ты мень­ше еди­ни­цы, то к кон­цу ра­бо­ты ал­го­рит­ма чис­ла ста­нут очень ма­лы.
  • По­лу­чен­ная тре­уголь­ная мат­ри­ца не при­год­на для про­сто­го вы­чис­ле­ния опре­де­ли­те­ля ис­ход­ной мат­ри­цы. Ибо, как не­слож­но по­счи­тать, опре­де­ли­тель D_T тре­уголь­ной мат­ри­цы те­перь ра­вен опре­де­ли­те­лю D_R ис­ход­ной мат­ри­цы, умно­жен­но­му на \left(a_{00}\right)^{N-1}\cdot\left(a_{11}\right)^{N-2}\cdot...\cdot\left(a_{N-2,N-2}\right)^1 — диа­го­наль­ные эле­мен­ты ре­зуль­ти­рую­щей мат­ри­цы.

Ба­рейс (Bareiss) до­ка­зал, что пра­вая часть вы­ра­же­ния (10) мо­жет быть по­де­ле­на на­це­ло на a_{j-1,j-1} (а не на a_{jj}, как в ал­го­рит­ме Гаус­са). При этом пред­по­ла­га­ет­ся, что a_{-1,-1}\equiv 1:

\vec a_i \leftarrow \frac{a_{jj}\cdot\vec a_i-a_{ij}\cdot\vec a_j}{a_{j-1,j-1}}.
(11)

Под «де­ле­ни­ем на це­ло» здесь по­ни­ма­ет­ся не то, что ре­зуль­тат де­ле­ния — це­лое чис­ло, а то, что не по­тре­бу­ет­ся бес­ко­неч­но­го чис­ла бит для хра­не­ния ре­зуль­та­та, и ре­зуль­тат де­ле­ния по­тре­бу­ет мень­ше бит, чем тре­бо­ва­лось для хра­не­ния чис­ли­те­ля. Ес­ли же ис­ход­ные эле­мен­ты мат­ри­цы — це­лые, то это озна­ча­ет, что ре­зуль­тат де­ле­ния — це­лое чис­ло.

Ещё од­ной ин­те­рес­ной осо­бен­но­стью ал­го­рит­ма яв­ля­ет­ся то, что в нём все­го од­но «не­ском­пен­си­ро­ван­ное» умно­же­ние для каж­до­го диа­го­наль­но­го эле­мен­та, то есть опре­де­ли­тель D_T ре­зуль­ти­рую­щей тре­уголь­ной мат­ри­цы ра­вен опре­де­ли­те­лю D_R ис­ход­ной мат­ри­цы, умно­жен­но­му на a_{00}\cdot a_{11}\cdot ...\cdot a_{N-2,N-2} (здесь име­ют­ся в ви­ду диа­го­наль­ные эле­мен­ты ре­зуль­ти­рую­щей мат­ри­цы):

D_T=D_R\times a_{00}\cdot a_{11}\cdot ...\cdot a_{N-2,N-2},
(12)

но

D_T=a_{00}\cdot a_{11}\cdot ...\cdot a_{N-1,N-1},
(13)

по­это­му

D_R=\frac{a_{00}\cdot a_{11}\cdot ...\cdot a_{N-1,N-1}}{a_{00}\cdot a_{11}\cdot ...\cdot a_{N-2,N-2}}=a_{N-1,N-1}.
(14)

То есть опре­де­ли­тель ис­ход­ной мат­ри­цы ра­вен пра­во­му ниж­не­му эле­мен­ту мат­ри­цы, по­лу­чен­ной в ре­зуль­та­те ал­го­рит­ма Ба­рей­са!

Фор­маль­ная за­пись ал­го­рит­ма Ба­рей­са:

\begin{array}{l}
a_{-1,-1}\equiv 1\\
\text{For }j=0,...,N-2\\
\left|\quad\begin{array}{l}
	\text{For }i=j+1,...,N-1\\
	\left|\quad\vec a_i \leftarrow \cfrac{a_{jj}\cdot\vec a_i-a_{ij}\cdot\vec a_j}{a_{j-1,j-1}}\right.
\end{array}\right.\\
\end{array}
(15)

Ал­го­ритм Ба­рей­са так же, как и ал­го­ритм Гаус­са, мо­жет быть улуч­шен пу­тём до­бав­ле­ния пе­ре­ста­нов­ки строк.

При­мер ра­бо­ты ал­го­рит­ма Ба­рей­са

Рас­смот­рим ра­бо­ту ал­го­рит­ма на при­ме­ре мат­ри­цы (3). Умно­жа­ем первую, вто­рую, и тре­тью стро́ки на 3, и от­ни­ма­ем от них ну­ле­вую стро­ку, умно­жен­ную на 4, 1, и 2 со­от­вет­ствен­но. Де­ле­ние на этом ша­ге не нуж­но, так как a_{-1,-1}\equiv 1. По­лу­чим:

\begin{pmatrix}3&2&3&4\\
0&4&-3&-10\\
0&10&9&5\\
0&5&-3&-5\end{pmatrix}
(16)

Те­перь вы­чи­та­ем из вто­рой и тре­тьей строк, умно­жен­ных на 4, первую стро­ку, умно­жен­ную на 10 и 5 со­от­вет­ствен­но (стро́ки ну­ме­ру­ют­ся с ну­ля). По­сле это­го де­лим стро́ки на 3:

\begin{pmatrix}3&2&3&4\\
0&4&-3&-10\\
0&0&22&40\\
0&0&1&10\end{pmatrix}
(17)

На­ко­нец, умно­жа­ем по­след­нюю стро­ку на 22, вы­чи­та­ем из неё пред­по­след­нюю, и де­лим ре­зуль­тат на 4:

\begin{pmatrix}3&2&3&4\\
0&4&-3&-10\\
0&0&22&40\\
0&0&0&45\end{pmatrix}
(18)

Об­ра­ти­те вни­ма­ние на ниж­ний пра­вый эле­мент: это опре­де­ли­тель ис­ход­ной мат­ри­цы. И нам да­же не по­тре­бо­ва­лись дро­би или чис­ла с пла­ваю­щей точ­кой для при­ве­де­ния мат­ри­цы к тре­уголь­но­му ви­ду и вы­чис­ле­ния это­го опре­де­ли­те­ля.

Про­грамм­ная реа­ли­за­ция ал­го­рит­ма Ба­рей­са

Реа­ли­за­ция ал­го­рит­ма Ба­рей­са на язы­ке про­грам­ми­ро­ва­ния C++ по­чти не от­ли­ча­ет­ся от со­от­вет­ствую­щей реа­ли­за­ции ал­го­рит­ма Гаус­са. От­ли­чия по­ка­за­ны вос­кли­ца­тель­ны­ми зна­ка­ми:

template<class T> T BareissDeterminant(T **const M, int const N)
{
    T denom(1); //!
    int exchanges(0);
   
    for(int l1(0); l1<N-1; ++l1)
    { //Перебираю все строки матрицы, кроме последней
        int maxN( l1 );
        T maxValue( fabs( M[l1][l1] ) );
        for(int l2(l1+1); l2<N; ++l2)
        { //Нахожу строку с максимальным по модулю элементом
            T const value( fabs( M[l2][l1] ) );
            if( value > maxValue ) { maxN=l2; maxValue=value; }
        }
       
        if( maxN > l1 )
        { //Нужен обмен
            T *const temp(M[l1]); M[l1]=M[maxN]; M[maxN]=temp;
            ++exchanges;
        } else { //Обмен не нужен, но нужна проверка на ноль
            if( maxValue == T(0) ) return maxValue;
        }

        T const value1( M[l1][l1] ); //!
          //!
       
        for(int l2(l1+1); l2<N; ++l2)
        { //Вычитаю строку из всех последующих
            T const value2( M[l2][l1] ); //!
            M[l2][l1] = T(0);
            for(int c(l1+1); c<N; ++c) //!
                M[l2][c]=(M[l2][c]*value1-M[l1][c]*value2)/denom;
        }

        denom = value1; //!
    }
    //!
    if(exchanges%2) return -M[N-1][N-1]; else return M[N-1][N-1];
}

Ли­стинг 2. Реа­ли­за­ция ал­го­рит­ма Ба­рей­са с пе­ре­ста­нов­кой строк

Про­вер­ка ал­го­рит­мов на мат­ри­це Гиль­бер­та

Мат­ри­цей Гиль­бер­та на­зы­ва­ет­ся мат­ри­ца ви­да:

H_{ij}=\frac{1}{i+j+1},\quad i,j=0,...,N-1;
(19)

или

H^N=\begin{pmatrix}1&1/2&\ldots &1/N\\
1/2&1/3&\ldots &1/(N+1)\\
\vdots &\vdots &\ddots &\vdots\\
1/N&1/(N+1)&\ldots &1/(2N-1)\end{pmatrix}.
(20)

Мат­ри­ца Гиль­бер­та — клас­си­че­ский при­мер пло­хо обу­слов­лен­ной мат­ри­цы. Си­сте­мы урав­не­ний с мат­ри­цей Гиль­бер­та пло­хо ре­ша­ют­ся чис­лен­ны­ми ме­то­да­ми. Мы же по­про­бу­ем по­счи­тать опре­де­ли­тель мат­ри­цы Гиль­бер­та де­ся­то­го по­ряд­ка.

Опре­де­ли­те­ли об­рат­ной мат­ри­цы Гиль­бер­та об­ра­зу­ют сле­ду­ю­щую по­сле­до­ва­тель­ность чи­сел (опре­де­ли­тель об­рат­ной мат­ри­цы ра­вен об­рат­ной ве­ли­чи­не опре­де­ли­те­ля):

\begin{aligned}
\left|H^1\right|^{-1}&=1,\\
\left|H^2\right|^{-1}&=12,\\
\left|H^3\right|^{-1}&=2160,\\
\left|H^4\right|^{-1}&=6048000,\\
\left|H^5\right|^{-1}&=266716800000,\\
\left|H^6\right|^{-1}&=186313420339200000,\\
\left|H^7\right|^{-1}&=2067909047925770649600000,\\
\left|H^8\right|^{-1}&=365356847125734485878112256000000,\\
\left|H^9\right|^{-1}&=1028781784378569697887052962909388800000000,\\
\left|H^{10}\right|^{-1}&={}\\
{}=46206&89394791469131629562883903627872698368000000000.
\end{aligned}
(21)
\begin{aligned}
\left|H^1\right|^{-1}&=1,\\
\left|H^2\right|^{-1}&=12,\\
\left|H^3\right|^{-1}&=2160,\\
\left|H^4\right|^{-1}&=6048000,\\
\left|H^5\right|^{-1}&=266716800000,\\
\left|H^6\right|^{-1}&=186313420339200000,\\
\left|H^7\right|^{-1}&=2067909047925770649600000,\\
\left|H^8\right|^{-1}&=365356847125734485878112256000000,\\
\left|H^9\right|^{-1}&=1028781784378569697887052962909388800000000,\\
\left|H^{10}\right|^{-1}&=4620689394791469131629562883903627872698368000000000.
\end{aligned}
(21)

По­про­бу­ем най­ти \left|H^{10}\right|^{-1} с по­мо­щью опи­сан­ных ал­го­рит­мов, ис­поль­зуя чис­ла двой­ной точ­но­сти:

Без пе­ре­ста­нов­ки строкС пе­ре­ста­нов­кой строк
Ал­го­ритм Гаус­са4.620287...\cdot 10^{52}4.620207...\cdot 10^{52}
Ал­го­ритм Ба­рей­са4.620204...\cdot 10^{52}4.620375...\cdot 10^{52}

Таб­ли­ца 1. Об­рат­ная ве­ли­чи­на опре­де­ли­те­ля мат­ри­цы Гиль­бер­та,
вы­чис­лен­но­го раз­лич­ны­ми спо­со­ба­ми

Не­смот­ря на то, что ал­го­ритм Ба­рей­са поз­во­ля­ет точ­но и быст­ро по­счи­тать опре­де­ли­тель мат­ри­цы, он ма­ло при­го­ден для ре­ше­ния си­стем ли­ней­ных урав­не­ний боль­шой раз­мер­но­сти. Де­ло в том, что по ме­ре при­ве­де­ния мат­ри­цы к тре­уголь­но­му ви­ду ме­то­дом Ба­рей­са, по­лу­чае­мые чис­ла по по­ряд­ку ве­ли­чи­ны рав­ны опре­де­ли­те­лю об­ра­бо­тан­ной на дан­ный мо­мент ча­сти мат­ри­цы. Мо­жет так ока­зать­ся, что опре­де­ли­те­ли на­столь­ко ве­ли­ки, что не уме­ща­ют­ся в ис­поль­зуе­мый тип дан­ных. Это озна­ча­ет, что опре­де­ли­тель та­кой мат­ри­цы не мо­жет быть по­счи­тан с ис­поль­зо­ва­ни­ем имею­ще­го­ся ти­па дан­ных, но это со­вер­шен­но не озна­ча­ет, что та­кую си­сте­му урав­не­ний нель­зя ре­шить. В та­ких слу­ча­ях луч­ше ис­поль­зо­вать ал­го­ритм Гаус­са с пе­ре­ста­нов­кой строк.

Мы ви­дим, что точ­нее все­го опре­де­ли­тель вы­чис­лен при по­мо­щи ал­го­ритм Ба­рей­са с пе­ре­ста­нов­ка­ми строк, хо­тя и ему да­ле­ко до пра­виль­но­го от­ве­та. Воз­мож­но, од­ной из при­чин столь зна­чи­тель­но­го от­кло­не­ния от пра­виль­но­го зна­че­ния яв­ля­ет­ся тот факт, что са­ми эле­мен­ты мат­ри­цы Гиль­бер­та бы­ли округ­ле­ны до чи­сел двой­ной точ­но­сти при ге­не­ра­ции мат­ри­цы.

Па­ра­л­лель­ная реа­ли­за­ция

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

Ес­ли посмот­реть ли­стинг 2, то мож­но за­ме­тить, что «вы­чи­та­ние» строки́ из всех по­сле­дую­щих — не­за­ви­си­мые опе­ра­ции, и об­ра­бот­ка эле­мен­тов каж­дой строки́ — то­же не­за­ви­си­мые опе­ра­ции. Но име­ет­ся за­ви­си­мость об­ра­бот­ки эле­мен­тов од­ной строки́ от об­ра­бот­ки эле­мен­тов дру­гой строки́ че­рез пе­ре­мен­ную value2.

Мож­но ли­бо да­вать каж­до­му вы­чис­ли­те­лю на об­ра­бот­ку це­лое чис­ло строк и свой эк­зем­пляр пе­ре­мен­ной value2 (что мо­жет при­ве­сти к дис­ба­лан­су на­груз­ки, так как чис­ло строк мо­жет не по­де­лить­ся на­це­ло на чис­ло вы­чис­ли­те­лей), ли­бо во­об­ще не ис­поль­зо­вать пе­ре­мен­ную value2 (а «пер­вые» эле­мен­ты строк об­ну­лять в кон­це, или во­об­ще не об­ну­лять, — на зна­че­ние опре­де­ли­те­ля это не по­вли­я­ет), и раз­бить всё мно­же­ство опе­ра­ций «вы­чи­та­ния» по­ров­ну меж­ду про­цес­со­ра­ми. В по­след­нем слу­чае дис­ба­ланс на­груз­ки ес­ли и воз­ник­нет, то это бу­дет не од­на стро­ка, как при раз­би­е­нии по стро­кам, а один эле­мент, что го­раз­до ме­нее за­мет­но. За­ве­до­мо ну­ле­вые эле­мен­ты век­то­ров не об­ра­ба­ты­ва­ют­ся, по­это­му об­щее чис­ло опе­ра­ций на этом эта­пе рав­но 4\cdot \left(N-l_1-1\right)^2 штук.

А вот на­хож­де­ние строки́ с мак­си­маль­ным по мо­ду­лю эле­мен­том — за­ви­си­мые опе­ра­ции (они опе­ри­ру­ют об­щи­ми пе­ре­мен­ны­ми maxN и maxValue). Мож­но, ко­неч­но, при­ме­нить в этом слу­чае ал­го­ритм сдваи­ва­ния, но го­раз­до про­ще оста­вить эту часть ал­го­рит­ма по­сле­до­ва­тель­ной, ибо на­хож­де­ние мак­си­маль­но­го эле­мен­та — ма­лая часть опе­ра­ций (при­мер­но 4\cdot \left(N-l_1-1\right) штук) по срав­не­нию с вы­чи­та­ни­ем строк.

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

Най­дём па­ра­метр \alpha для за­ко­на Ам­да­ля. По опре­де­ле­нию, \alpha ра­вен от­но­ше­нию ко­ли­че­ства по­сле­до­ва­тель­ных опе­ра­ций к об­ще­му чис­лу опе­ра­ций:

\alpha=\frac{N_1}{N_1+N_p}.
(22)

Чис­ло по­сле­до­ва­тель­ных опе­ра­ций:

N_1 = \sum_{l_1=0}^{N-2} 4\cdot \left(N-l_1-1\right) = 4\cdot\frac{N^2-N}{2}.
(23)

Чис­ло па­ра­л­лель­ных опе­ра­ций:

N_p = \sum_{l_1=0}^{N-2} 4\cdot \left(N-l_1-1\right)^2 = 4\cdot\frac{2N^3-3N^2+N}{6}.
(24)

Под­став­ляя вы­ра­же­ния (23) и (24) в (22), по­лу­чим:

\alpha\left(N\right)=\frac{3}{2N+2}.
(25)

По­лу­чен­ное вы­ра­же­ние мож­но под­ста­вить в вы­ра­же­ние (7) преды­ду­щей лек­ции, и ис­поль­зо­вать его, как оцен­ку свер­ху для воз­мож­но­го уско­ре­ния ал­го­рит­ма. Вид­но, что чем боль­ше N, тем мень­ше бу­дет \alpha, и тем бли­же кри­вая уско­ре­ния бу­дет к пря­мой мак­си­маль­но­го уско­ре­ния S_p=p. Эта за­ко­но­мер­ность ча­сто встре­ча­ет­ся на прак­ти­ке: чем боль­ше ис­ход­ных дан­ных да­ёт­ся ал­го­рит­му, тем обыч­но с боль­шей эф­фек­тив­но­стью этот ал­го­ритм мож­но рас­па­рал­ле­лить.

Что по­чи­тать по те­ме

Chee K. Yap, “Fundamental Problems in Algorithmic Algebra — лек­ции по чис­лен­ным ал­го­рит­мам ал­геб­ры, с до­ка­затель­ства­ми всех утвер­жде­ний. Ал­го­ритм Ба­рей­са опи­сан в лек­ции 10 (.ps-фай­лы бо­лее но­вые — чи­тай­те их).

27 отзывов на «приведение матрицы к треугольному виду»

От­лич­ный кон­спект. Един­ствен­ное, еще бы ме­тод Ха­лец­ко­го сю­да до­ба­вить — бы­ло бы во­об­ще СУПЕР!!!!!!!!!
а я так ни че­го и не по­нял!=(
Здрав­ствуй­те, ес­ли я ре­шаю си­сте­му ли­ней­ных од­но­род­ных урав­не­ний, то то­гда мне на­до, что­бы опре­де­ли­тель как раз рав­нял­ся 0. Мож­но ли с по­мо­щью это­го ме­то­да опре­де­лить, ра­вен ли опре­де­ли­тель 0 (с не­ко­то­рой сте­пе­нью точ­но­сти)?
На са­мом де­ле это очень слож­ный во­прос. Что́ счи­тать ра­вен­ством или не ра­вен­ством, ес­ли вы­чис­ле­ния вы­пол­ня­ют­ся при­бли­жён­но? Ес­ли у вас все ко­эф­фи­ци­ен­ты в ис­ход­ной мат­ри­це — це­лые чис­ла, то вам очень по­вез­ло — бе­ри­те ал­го­ритм Ба­рей­са с пе­ре­ста­нов­кой строк (ли­стинг 2) и, ес­ли функ­ция вер­ну­ла 0, то у вас опре­де­ли­тель в точ­но­сти ра­вен ну­лю. Ес­ли ис­ход­ные чис­ла ра­цио­наль­ные — то вы­пол­няй­те все вы­чис­ле­ния в дро­бях, и полу́чите точ­ный ре­зуль­тат. Ес­ли же ва­ши чис­ла — про­из­воль­ные чис­ла с пла­ваю­щей точ­кой, то си­ту­а­ция го­раз­до ху­же. Ведь вам для ре­ше­ния нуж­но при­ве­сти мат­ри­цу не к тре­уголь­ной, а к тра­пе­це­вид­ной фор­ме: то есть при об­на­ру­же­нии ну­ле­во­го под­столб­ца нуж­но сме­щать­ся впра­во на один стол­бец. Каж­дое та­кое «сме­ще­ние впра­во» умень­ша­ет ранг мат­ри­цы на еди­ни­цу. Ранг — это це­лое чис­ло, ко­то­рое вам нуж­но знать точ­но, ибо ранг опре­де­ля­ет раз­мер­ность под­про­стран­ства ре­ше­ний СЛАУ. Ма­лей­шие по­греш­но­сти при вы­чис­ле­ни­ях мо­гут при­ве­сти к ка­ко­му угод­но из­ме­не­нию ран­га мат­ри­цы. По­это­му со­вер­шен­но не­по­нят­но, что де­лать в этом слу­чае. Что́ счи­тать чис­лом, близ­ким к ну­лю? Воз­мож­но, сле­ду­ет по­сту­пить сле­дую­щим об­ра­зом: най­ти в ис­ход­ной мат­ри­це мак­си­маль­ное по мо­ду­лю чис­ло, и мо­дуль это­го чис­ла при­нять за эта­лон. За­тем вос­поль­зо­вать­ся ме­то­дом Гаус­са с пе­ре­ста­нов­кой строк (этот ме­тод в сред­нем не ме­ня­ет ве­ли­чи­ну чи­сел в стро­ках мат­ри­цы). И, ес­ли в ка­кой-то мо­мент те­ку­щий зна­ме­натель ока­зы­ва­ет­ся по мо­ду­лю в не­сколь­ко мил­ли­о­нов раз мень­ше эта­ло­на, то он счи­та­ет­ся ну­лём, и ал­го­ритм пе­ре­хо­дит на один стол­бец пра­вее, тем са­мым умень­шая ранг на еди­ни­цу.
А ка­кой(или чей) ме­тод на­хож­де­ния опре­де­ли­те­ля мат­ри­цы са­мый точ­ный?
Всё за­ви­сит от то­го, в ка­ком ви­де у вас пред­став­ле­ны дан­ные. Ес­ли дан­ные у вас пред­став­ле­ны в це­лых чис­лах, и вы­чис­ля­е­те вы в це­лых чис­лах с до­ста­точ­ным диа­па­зо­ном, то ал­го­ритм Ба­рей­са даст точ­ный от­вет (не важ­но, с пе­ре­ста­нов­кой строк или без, лишь бы де­ле­ния на ноль не бы­ло). Ес­ли дан­ные пред­став­ле­ны в дро­бях, и вы­чис­ле­ния вы ве­дё­те в дро­бях, то все ме­то­ды да­дут точ­ный от­вет. Ес­ли же у вас чис­ла с пла­ваю­щей точ­кой огра­ни­чен­ной точ­но­сти (и ис­ход­ные дан­ные то­го же ти­па), то са­мый точ­ный от­вет бу­дет у ал­го­рит­ма Ба­рей­са с пе­ре­ста­нов­кой строк.
ОГРОМНЕШОЕ СПАСИБО, осо­бен­но за реа­ли­за­цию ал­го­рит­ма Гаус­са
3 дня про­му­чал­ся над ним а ощиб­ка бы­ла все­го в не­до­стою­щей 1-й стро­ки про­ги ещё раз спа­си­бо
Спа­си­бо боль­шое за лек­цию, очень по­мог­ла при ре­ше­нии прак­ти­ки сту­ден­ту-ма­те­ма­ти­ку)
Чув­ствую се­бя чай­ни­ком, я и про­стые мат­ри­цы не мо­гу при­ве­сти к тре­уг. ви­ду(((( пре­под­ша ме­ня зав­тра съест без со­ли(((((
Здрав­ствуй­те! Спа­си­бо за столь хо­ро­ший кон­спект! Очень по­мог­ло. Но не до кон­ца( Учусь про­грам­ми­ро­вать в СИ, на­до на­пи­сать про­грам­му, ко­то­рая бу­дет ре­шать СЛАУ ме­то­дом Гаус­са. В прин­ци­пе по­чти сде­лал, но ни­как не мо­гу на­пи­сать усло­вие, ко­гда мат­ри­ца не име­ет ре­ше­ния или ко­гда бес­ко­неч­ное мно­же­ство ре­ше­ний. Ес­ли кто мо­жет — по­мо­ги­те по­жа­луй­ста. Ски­ды­ваю код в СИ:
#include
#include
#define N 20
#define M (N+1) void glavelem( int k, double mas[][N + 1], int n, int otv[] )
{
  int i, j, i_max = k, j_max = k;
  double temp;
  //Ищем максимальный по модулю элемент
  for ( i = k; i < n; i++ )
    for ( j = k; j < n; j++ )
      if ( fabs( mas[i_max] [j_max] ) < fabs( mas[i] [j] ) )
      {
        i_max = i;
        j_max = j;
      }
  //Переставляем строки
  for ( j = k; j < n + 1; j++ )
  {
    temp = mas[k] [j];
    mas[k] [j] = mas[i_max] [j];
    mas[i_max] [j] = temp;
  }
  //Переставляем столбцы
  for ( i = 0; i N )
    {
      printf("Слишком большое кол-во переменных.n");
      return 0;
    }
  }
  while ( N < n );   printf( "Задайте %d элементов матрицы:n", n*(n+1));
  for ( i = 0; i < n; i++ )
    for ( j = 0; j < n + 1; j++ )
      scanf( "%lf", &amp;mas[i][j] );
     
  /*——————————————————*/
  //Сначала все корни по порядку
  for ( i = 0; i < n + 1; i++ )
    otv[i] = i;
  //Прямой ход метода Гаусса
  for ( k = 0; k = k; j )
      mas[k] [j] /= mas[k] [k];
    for ( i = k + 1; i = k; j )
        mas[i] [j] -= mas[k] [j] * mas[i] [k];
  }  
  //Обратный ход
  for ( i = 0; i = 0; i )
    for ( j = i + 1; j < n; j++ )
      x[i] -= x[j] * mas[i] [j];
  //Вывод результата
  printf( "Одно решение:n" );
  for ( i = 0; i < n; i++ )
    for ( j = 0; j < n; j++ )
      if ( i == otv[j] )
      { //Расставляем корни по порядку
        printf( "x%d = %fn", y++, x[j] );
        break;
      }
  return ( 0 );
}
П.С за­ра­нее спа­си­бо тем, кто от­клик­нет­ся. Ра­бо­ту сда­вать зав­тра. По воз­мож­но­сти на­пи­ши­те не­до­стаю­щие стро­ки ко­да:
1)ес­ли мат­ри­ца не име­ет ре­ше­ния
2)ес­ли мат­ри­ца име­ет бес­ко­неч­ное мно­же­ство ре­ше­ний
не знаю по­че­му, но часть ко­да про­гла­ты­ва­ет… но ес­ли вдруг кто то за­хо­чет по­мочь — пи­ши­те на ящик: apelsin_ast@mail.ru
На­де­юсь на Ва­шу по­мощь
<< "(прак­ти­че­ской реа­ли­за­ци­ей мы зай­мём­ся в од­ной из сле­дую­щих лек­ций)." По­че­му нет обе­щан­ной па­ра­л­лель­ной реа­ли­за­ции при­ве­де­ния к тре­уголь­но­му ви­ду?
Это наи­бо­лее ин­те­ре­су­ет.
Спа­си­бо, от­лич­ный ма­те­ри­ал..очень по­мог­ли!
от­лич­ная лек­ция, спа­си­бо)
остал­ся во­прос….Оба ме­то­да мож­но на эк­за­ме­не ис­поль­зо­вать?или ка­кой-то один?
В функ­ции BareissDeterminant не учли слу­чай ко­гда N = 1, то есть ко­гда мат­ри­ца име­ет все­го один эле­мент.
Нет, в ко­де всё вер­но. При N = 1 от­ра­бо­та­ет return M[N-1][N-1];
А, разо­брал­ся, спа­си­бо!
Гор­но­лыж­ные ко­стю­мы и спор­тив­ную брен­до­вую одеж­ду, обувь и ак­сес­су­а­ры для муж­чин, жен­щин и де­тей по са­мым низ­ким це­нам ку­пить спор­тив­ные ко­стю­мы, Зим­ние бейс­бол­ки и мно­гое дру­гое в мар­кет-мар­кет.ру http://market-market.ru
гру­зо­пе­ре­воз­ки в ниж­нем нов­го­ро­де http://perevozkinn.com/
ин­ди­ви­ду­ал­ки но­во­си­бир­ска http://pikapnsk.com/
Про­сти­тут­ки НИЖНЕГО Нов­го­ро­да http://privatnnover.com лю­бов­ни­цы и ин­ди­ви­ду­ал­ки. Зна­чи­мые со­бы­тия для зна­ком­ства и раз­вле­че­ний.
ба­ни ниж­не­го нов­го­ро­да сау­на Эдем Ниж­ний Нов­го­род
В се­ти об­на­ру­же­на но­вая си­сте­ма по за­ра­бот­ку он­лайн. http://magstreet.ru/ — ра­бо­та без вло­же­ний Га­ран­ти­ро­ван­ная при­быль.Соб­ствен­ный он­лайн сер­вис
ноч­ная до­став­ка ал­ко­го­ля до­став­ка вод­ки но­чью Се­ва­сто­поль
ин­ди­ви­ду­ал­ки ниж­не­го нов­го­ро­да про­сти­тут­ки ниж­ний нов­го­род гос­по­жа
Я во­об­ще ни­че­го не по­ни­маю в Ва­шей та­ра­бар­щи­не, я гу­ма­ни­та­рист, ну Вы и мо­лод­цы, как мож­но эту всю фиг­ню по­нять , да еще и при­ме­нить ку­да-ни­будь. Я вос­хи­ще­на про­сто. По­мо­гаю пи­сать ди­плом на те­му Гео­мет­рия тре­уголь­ной мат­ри­цы. Вот ищу в ин­тер­не­те, для ме­ня это но­во­сти, до мо­е­го ума не со­всем до­хо­дя­щие. Ну, мо­лод­цы!!!
Rattling great visual appeal on this web site, I’d value it 10.

Оставить отзыв

Жёлтые поля обязательны к заполнению

   

Можете использовать теги <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> <pre lang=""> <div class=""> <span class=""> <br>