chap31.htm

介绍算法的书籍

<h2><a name="0979_1ab0">Exercises<a name="0979_1ab0"></h2><P>
<a name="0979_1ab1">31.3-1<a name="0979_1ab1"><P>
Does Strassen's algorithm work over the number system (<B>Z</B>[<I>x</I>], +, <img src="../images/dot10.gif">, 0, 1), where <B>Z</B>[<I>x</I>] is the set of all polynomials with integer coefficients in the variable <I>x</I> and + and <img src="../images/dot10.gif"> are ordinary polynomial addition and multiplication?<P>
<a name="0979_1ab2">31.3-2<a name="0979_1ab2"><P>
Explain why Strassen's algorithm doesn't work over closed semirings (see Section 26.4) or over the boolean quasiring ({0, 1}, <FONT FACE="Times New Roman" SIZE=4><img src="../images/angledwn.gif"></FONT>, <FONT FACE="Times New Roman" SIZE=4><img src="../images/angleup.gif"></FONT>, 0, 1).<P>
<a name="0979_1ab3">31.3-3<a name="0979_1ab3"><P>
Prove Theorem 31.7 and Corollary 31.8.<P>
<a name="0979_1ab4">31.3-4<a name="0979_1ab4"><P>
Show that the boolean quasiring ({0, 1}, <FONT FACE="Times New Roman" SIZE=4><img src="../images/angledwn.gif"></FONT>, <FONT FACE="Times New Roman" SIZE=4><img src="../images/angleup.gif"></FONT>, 0, 1) cannot be embedded in a ring. That is, show that it is impossible to add a "-1" to the quasiring so that the resulting algebraic structure is a ring.<P>
<a name="0979_1ab5">31.3-5<a name="0979_1ab5"><P>
Argue that if all computations in the algorithm of Theorem 31.10 are performed modulo <I>n</I> + 1, the algorithm still works correctly.<P>
<a name="0979_1ab6">31.3-6<a name="0979_1ab6"><P>
<a name="0979_1aab">Show how to determine efficiently if a given undirected input graph contains a triangle (a set of three mutually adjacent vertices).<P>
<a name="0979_1ab7">31.3-7<a name="0979_1ab7"><P>
<a name="0979_1aac"><a name="0979_1aad"><a name="0979_1aae"><a name="0979_1aaf">Show that computing the product of two <I>n</I> <img src="../images/mult.gif"> <I>n</I> boolean matrices over the boolean quasiring is reducible to computing the transitive closure of a given directed 3<I>n</I>-vertex input graph.<P>
<a name="0979_1ab8">31.3-8<a name="0979_1ab8"><P>
Show how to compute the transitive closure of a given directed <I>n</I>-vertex input graph in time <I>O</I>(<I>n</I><SUP>lg 7</SUP> lg <I>n</I>). Compare this result with the performance of the <FONT FACE="Courier New" SIZE=2>TRANSITIVE</FONT>-<FONT FACE="Courier New" SIZE=2>CLOSURE</FONT> procedure in Section 26.2.<P>
<P>


<P>







<h1><a name="097a_1ab5">31.4 Solving systems of linear equations<a name="097a_1ab5"></h1><P>
<a name="097a_1ab0"><a name="097a_1ab1">Solving a set of simultaneous linear equations is a fundamental problem that occurs in diverse applications. A linear system can be expressed as a matrix equation in which each matrix or vector element belongs to a field, typically the real numbers <B>R</B>. This section discusses how to solve a system of linear equations using a method called LUP decomposition.<P>
We start with a set of linear equations in <I>n</I> unknowns <I>x</I><SUB>1</SUB>,<I> x</I><SUB>2</SUB>, . . . ,<I> x<SUB>n</I></SUB>:<P>
<pre><I>a</I><SUB>11</SUB><I>x</I><SUB>1</SUB> + <I>a</I><SUB>12</SUB><I>x</I><SUB>2</SUB> + <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> + <I>a</I><SUB>1<I>n</SUB>x<SUB>n</I></SUB> = <I>b</I><SUB>1</SUB>,</sub></sup></pre><P>
<pre><I>a</I><SUB>21</SUB><I>x</I><SUB>1</SUB> + <I>a</I><SUB>22</SUB><I>x</I><SUB>2</SUB> + <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> + <I>a</I><SUB>2<I>n</SUB>x<SUB>n</SUB> </I>= <I>b</I><SUB>2</SUB>,</sub></sup></pre><P>
<img src="750_a.gif"><P>
<h4><a name="097a_1ab6">(31.17)<a name="097a_1ab6"></sub></sup></h4><P>
<pre><I>a<SUB>n</I>1</SUB><I>x</I><SUB>1</SUB> + <I>a<SUB>n</I>2</SUB><I>x</I><SUB>2</SUB> + <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> + <I>a<SUB>nn</SUB>x<SUB>n</SUB> </I>= <I>b<SUB>n</I></SUB>.</sub></sup></pre><P>
<a name="097a_1ab2">A set of values for <I>x</I><SUB>1</SUB>, <I>x</I><SUB>2</SUB>, . . . , <I>x<SUB>n</I></SUB> that satisfy all of the equations (31.17) simultaneously is said to be a <I><B>solution</I></B> to these equations. In this section, we only treat the case in which there are exactly <I>n</I> equations in <I>n</I> unknowns.<P>
We can conveniently rewrite equations (31.17) as the matrix-vector equation<P>
<img src="750_b.gif"><P>
or, equivalently, letting <I>A</I> = (<I>a<SUB>ij</I></SUB>), <I>x</I> = (<I>x<SUB>j</I></SUB>), and <I>b</I> = (<I>b<SUB>i</I></SUB>), as<P>
<pre><I>Ax</I> = <I>b </I>.</sub></sup></pre><P>
<h4><a name="097a_1ab7">(31.18)<a name="097a_1ab7"></sub></sup></h4><P>
If <I>A</I> is nonsingular, it possesses an inverse <I>A</I><SUP>-1</SUP>, and<P>
<pre><I>x</I> = <I>A</I><SUP>-1 </SUP><I>b</I></sub></sup></pre><P>
<h4><a name="097a_1ab8">(31.19)<a name="097a_1ab8"></sub></sup></h4><P>
is the solution vector. We can prove that <I>x</I> is the unique solution to equation (31.18) as follows. If there are two solutions, <I>x </I>and <I>x</I>', then <I>Ax</I> = <I>Ax</I>' = <I>b</I> and<P>
<pre><I>x  </I>=  (<I>A</I><SUP>-1 </SUP><I>A</I>)<I>x</I></sub></sup></pre><P>
<pre>=  <I>A</I><SUP>-1</SUP>(<I>Ax</I>)</sub></sup></pre><P>
<pre>=  <I>A</I><SUP>-1</SUP>(<I>Ax</I>')</sub></sup></pre><P>
<pre>=  (<I>A</I><SUP>-1 </SUP><I>A</I>)<I>x</I>'</sub></sup></pre><P>
<pre>=  <I>x</I>'.</sub></sup></pre><P>
<a name="097a_1ab3"><a name="097a_1ab4">In this section, we shall be concerned predominantly with the case in which <I>A</I> is nonsingular or, equivalently (by Theorem 31.1), the rank of <I>A</I> is equal to the number <I>n</I> of unknowns. There are other possibilities, however, which merit a brief discussion. If the number of equations is less than the number <I>n</I> of unknowns--or, more generally, if the rank of <I>A</I> is less than <I>n</I>--then the system is <I><B>underdetermined</I></B>. An underdetermined system typically has infinitely many solutions (see Exercise 31.4-9), although it may have no solutions at all if the equations are inconsistent. If the number of equations exceeds the number <I>n</I> of unknowns, the system is <I><B>overdetermined</I></B>, and there may not exist any solutions. Finding good approximate solutions to overdetermined systems of linear equations is an important problem that is addressed in Section 31.6.<P>
Let us return to our problem of solving the system<I> Ax</I> = <I>b</I> of <I>n</I> equations in <I>n</I> unknowns. One approach is to compute <I>A</I><SUP>-1</SUP> and then multiply both sides by <I>A</I><SUP>-1</SUP>, yielding <I>A</I><SUP>-1</SUP><I>Ax </I>=<I> A</I><SUP>-1</SUP><I>b</I>, or <I>x </I>=<I> A</I><SUP>-1</SUP><I>b</I>. This approach suffers in practice from <I><B>numerical instability</I></B>: round-off errors tend to accumulate unduly when floating-point number representations are used instead of ideal real numbers. There is, fortunately, another approach--LUP decomposition--that is numerically stable and has the further advantage of being about a factor of 3 faster.<P>





<h2>Overview of LUP decomposition</h2><P>
<a name="097b_1ab5"><a name="097b_1ab6"><a name="097b_1ab7">The idea behind LUP decomposition is to find three <I>n</I> <img src="../images/mult.gif"><I> n</I> matrices <I>L, U</I>, and <I>P</I> such that<P>
<pre><I>PA</I> = <I>LU </I>,</sub></sup></pre><P>
<h4><a name="097b_1ab8">(31.20)<a name="097b_1ab8"></sub></sup></h4><P>
where<P>
<FONT FACE="Courier New" SIZE=2><img src="../images/dot12.gif"></FONT>     <I>L</I> is a unit lower-triangular matrix,<P>
<FONT FACE="Courier New" SIZE=2><img src="../images/dot12.gif"></FONT>     <I>U</I> is an upper-triangular matrix, and<P>
<FONT FACE="Courier New" SIZE=2><img src="../images/dot12.gif"></FONT>     <I>P</I> is a permutation matrix.<P>
We call matrices <I>L, U,</I> and <I>P</I> satisfying equation (31.20) an <I><B>LUP decomposition</I></B> of the matrix <I>A</I>. We shall show that every nonsingular matrix <I>A</I> possesses such a decomposition.<P>
The advantage of computing an LUP decomposition for the matrix <I>A</I> is that linear systems can be solved more readily when they are triangular, as is the case for both matrices <I>L</I> and <I>U</I>. Having found an LUP decomposition for <I>A</I>, we can solve the equation (31.18) <I>Ax </I>=<I> b</I> by solving only triangular linear systems, as follows. Multiplying both sides of <I>Ax </I>=<I> b</I> by <I>P</I> yields the equivalent equation <I>P Ax </I>=<I> Pb</I>, which by Exercise 31.1-2 amounts to permuting the equations (31.17). Using our decomposition (31.20), we obtain<P>
<pre><I>LUx</I> = <I>Pb</I> .</sub></sup></pre><P>
We can now solve this equation by solving two triangular linear systems. Let us define <I>y </I>=<I> Ux</I>, where <I>x</I> is the desired solution vector. First, we solve the lower-triangular system<P>
<pre><I>Ly = Pb</I></sub></sup></pre><P>
<h4><a name="097b_1ab9">(31.21)<a name="097b_1ab9"></sub></sup></h4><P>
for the unknown vector <I>y</I> by a method called "forward substitution." Having solved for <I>y,</I> we then solve the upper-triangular system<P>
<pre><I>Ux = y</I></sub></sup></pre><P>
<h4><a name="097b_1aba">(31.22)<a name="097b_1aba"></sub></sup></h4><P>
for the unknown <I>x</I> by a method called "back substitution." The vector <I>x</I> is our solution to <I>Ax </I>=<I> b</I>, since the permutation matrix <I>P</I> is invertible (Exercise 31.1 -2):<P>
<pre><I>Ax  =  P<SUP>-</I>1</SUP> <I>LUx</I></sub></sup></pre><P>
<pre>=  <I>P<SUP>-</I>1</SUP> <I>Ly</I></sub></sup></pre><P>
<pre>=  <I>P<SUP>-</I>1</SUP> <I>Pb</I></sub></sup></pre><P>
<pre>=  <I>b</I> .</sub></sup></pre><P>
Our next step is to show how forward and back substitution work and then attack the problem of computing the LUP decomposition itself.<P>
<P>







<h2>Forward and back substitution</h2><P>
<a name="097c_1ab8"><I><B>Forward substitution</I></B> can solve the lower-triangular system (31.21) in <img src="../images/bound.gif">(<I>n</I><SUP>2</SUP>) time, given<I> L, P, </I>and <I>b</I>. For convenience, we represent the permutation <I>P</I> compactly by an array <img src="../images/piuc.gif">[1 . . <I>n</I>]. For <I>i</I> = 1, 2, . . . , <I>n</I>, the entry <img src="../images/piuc.gif">[<I>i</I>] indicates that <I>P<SUB>i</I>,</SUB><img src="../images/piuc.gif">[<I>i</I>] = 1 and <I>P<SUB><FONT FACE="Courier New" SIZE=2>ij</I></FONT></SUB> = 0 for<I> j </I><img src="../images/noteq.gif"><I> </I><img src="../images/piuc.gif">[<I>i</I>]. Thus, <I>PA</I> has <I>a</I><SUB><FONT FACE="Courier New" SIZE=2></SUB><img src="../images/piuc.gif">[<I>i</I>],<I>j<SUB></SUB> </I></FONT>in row <I>i </I>and column <I>j</I>, and<I> Pb</I> has <I>b</I><img src="../images/piuc.gif">[<I>i</I>]<SUB> </SUB>as its <I>i</I>th element. Since <I>L</I> is unit lower-triangular, equation (31.21) can be rewritten as<P>
<pre><I>y</I><SUB>1</SUB>                             =  <I>b</I><SUB></SUB><img src="../images/piuc.gif"><SUB>[1],</sub></sup></pre><P>
<pre><I>l</I><SUB>21</SUB><I>y</I><SUB>1</SUB> + <I>y</I><SUB>2                        </SUB>=  <I>b</I><SUB></SUB><img src="../images/piuc.gif"><SUB>[2],</sub></sup></pre><P>
<pre><I>l</I><SUB>31</SUB><I>y</I><SUB>1</SUB> + <I>l</I><SUB>32</SUB><I>y</I><SUB>2</SUB> + <I>y</I><SUB>3                </SUB>=  <I>b</I><SUB></SUB><img src="../images/piuc.gif"><SUB>[3],</sub></sup></pre><P>
<img src="752_a.gif"><P>
<pre><I>l<SUB>n</I>1</SUB><I>y</I><SUB>1</SUB> + <I>l<SUB>n</I>2</SUB><I>y</I><SUB>2</SUB> + <I>l<SUB>n</I>3</SUB><I>y</I><SUB>3</SUB> + <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> <img src="../images/dot10.gif"> + <I>y<SUB>n  </I></SUB>=  <I>b</I><SUB></SUB><img src="../images/piuc.gif"><SUB>[<I>n</I>]</SUB>.</sub></sup></pre><P>
Quite apparently, we can solve for <I>y</I><SUB>1</SUB> directly, since the first equation tells us that <I>y</I><SUB>1</SUB> = <I>b</I><img src="../images/piuc.gif"><SUB>[1]</SUB>. Having solved for <I>y</I><SUB>1</SUB>, we can substitute it into the second equation, yielding<P>
<pre><I