Лек­ция № 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-фай­лы бо­лее но­вые — чи­тай­те их).

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

От­лич­ный кон­спект. Един­ствен­ное, еще бы ме­тод Ха­лец­ко­го сю­да до­ба­вить — бы­ло бы во­об­ще СУПЕР!!!!!!!!!
а я так ни че­го и не по­нял!=(
Здрав­ствуй­те, ес­ли я ре­шаю си­сте­му ли­ней­ных од­но­род­ных урав­не­ний, то то­гда мне на­до, что­бы опре­де­ли­тель как раз рав­нял­ся 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.
I’m not sure exactly why but this web site is loading extremely slow for me.
Is anyone else having this problem or is it a problem on my end?
I’ll check back later and see if the problem still exists.
If you wish for to increase your knowledge simply keep visiting this web page
and be updated with the most recent gossip posted here.
Thanks a lot for sharing this with all folks you really
recognize what you are speaking approximately! Bookmarked.
Please additionally consult with my website =). We can have a hyperlink exchange agreement between us

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

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

   

Можете использовать теги <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>