Linear Algebra/Topic: Speed of Calculating Determinants/Solutions

Solutions

Most of these problems presume access to a computer.

Problem 1

Computer systems generate random numbers (of course, these are only pseudo-random, in that they are generated by an algorithm, but they pass a number of reasonable statistical tests for randomness).

1. Fill a $5\!\times \!5$ array with random numbers (say, in the range $[0..1)$). See if it is singular. Repeat that experiment a few times. Are singular matrices frequent or rare (in this sense)?
2. Time your computer algebra system at finding the determinant of ten $5\!\times \!5$ arrays of random numbers. Find the average time per array. Repeat the prior item for $15\!\times \!15$ arrays, $25\!\times \!25$ arrays, and $35\!\times \!35$ arrays. (Notice that, when an array is singular, it can sometimes be found to be so quite quickly, for instance if the first row equals the second. In the light of your answer to the first part, do you expect that singular systems play a large role in your average?)
3. Graph the input size versus the average time.

Answer

1. Under Octave, rank(rand(5)) finds the rank of a $5\!\times \!5$ matrix whose entries are (uniformily distributed) in the interval $[0..1)$. This loop which runs the test $5000$ times octave:1> for i=1:5000
   > if rank(rand(5))<5 printf("That's one."); endif
   > endfor
   produces (after a few seconds) returns the prompt, with no output. The Octave script
   
   function elapsed_time = detspeed (size)
   a=rand(size);
   tic();
   for i=1:10
   det(a);
   endfor
   elapsed_time=toc();
   endfunction
   
   
   lead to this session.
   
   octave:1> detspeed(5)
   ans = 0.019505
   octave:2> detspeed(15)
   ans = 0.0054691
   octave:3> detspeed(25)
   ans = 0.0097431
   octave:4> detspeed(35)
   ans = 0.017398
   
2. Here is the data (rounded a bit), and the graph.

   ${\begin{array}{r|ccccccccc}{\textit {matrix\ rows}}&15&25&35&45&55&65&75&85&95\\\hline {\textit {time\ per\ ten}}&0.0034&0.0098&0.0675&0.0285&0.0443&0.0663&0.1428&0.2282&0.1686\end{array}}$

   (This data is from an average of twenty runs of the above script, because of the possibility that the randomly chosen matrix happens to take an unusually long or short time. Even so, the timing cannot be relied on too heavily; this is just an experiment.)

   

Problem 2

Compute the determinant of each of these by hand using the two methods discussed above.

1. ${\begin{vmatrix}2&1\\5&-3\end{vmatrix}}$
2. ${\begin{vmatrix}3&1&1\\-1&0&5\\-1&2&-2\end{vmatrix}}$
3. ${\begin{vmatrix}2&1&0&0\\1&3&2&0\\0&-1&-2&1\\0&0&-2&1\end{vmatrix}}$

Count the number of multiplications and divisions used in each case, for each of the methods. (On a computer, multiplications and divisions take much longer than additions and subtractions, so algorithm designers worry about them more.)

Answer

The number of operations depends on exactly how the operations are carried out.

1. The determinant is $-11$. To row reduce takes a single pivot with two multiplications ($-5/2$ times $2$ plus $5$, and $-5/2$ times $1$ plus $-3$) and the product down the diagonal takes one more multiplication. The permutation expansion takes two multiplications ($2$ times $-3$ and $5$ times $1$).
2. The determinant is $-39$. Counting the operations is routine.
3. The determinant is $4$.

Problem 3

What $10\!\times \!10$ array can you invent that takes your computer system the longest to reduce? The shortest?

Answer

One way to get started is to compare these under Octave: det(rand(10));, versus det(hilb(10));, versus det(eye(10));, versus det(zeroes(10));. You can time them as in tic(); det(rand(10)); toc().

Problem 4

Write the rest of the FORTRAN program to do a straightforward implementation of calculating determinants via Gauss' method. (Don't test for a zero pivot.) Compare the speed of your code to that used in your computer algebra system.

Answer

This is a simple one.

DO 5 ROW=1, N
PIVINV=1.0/A(ROW,ROW)
DO 10 I=ROW+1, N
DO 20 J=I, N
A(I,J)=A(I,J)-PIVINV*A(ROW,J)
20 CONTINUE
10 CONTINUE
5 CONTINUE

Problem 5

The FORTRAN language specification requires that arrays be stored "by column", that is, the entire first column is stored contiguously, then the second column, etc. Does the code fragment given take advantage of this, or can it be rewritten to make it faster, by taking advantage of the fact that computer fetches are faster from contiguous locations?

Answer

Yes, because the $J$ is in the innermost loop.