Programming Assignment #3: Repeated Squaring Solution

In section 3.7 of the text  we were introduced to the Power Method  for SVD computation. We had to deal with the problem of taking  an d × d matrix  B and raising it to a large exponent. That is, we have to compute  Bk for large values of k.

You could compute  it by multiplying B  over with  it  self k times.   That is, (with NEWMAT ) something  along the lines of

1:  C = B;

2:  for  int i = 1; i < k; i++ do

3:       C = B*C;

4:  end for

Assuming  you are doing straightforward matrix  multiplication (nothing clever like Strassen’s  method,  etc.)  then  the  above procedure  will take  O(d3 k) steps, as each matrix  multiplication is O(d3 ) and there  are k-many  of them.

A candidate algorithm that is more efficient than  the one shown above goes by the name of Repeated  Squaring,  and it described  below.

 

Repeated Squaring Algorithm

 

Suppose you want to compute  B11 , you write the exponent 11 in binary  – which is h1 0 1 1i2 . That is, 11 = 1 × 23 + 0 × 22 + 1 × 21 + 1 × 20 . You then  compute all dlog2 11e-many powers of B.  That is, you compute  B, B2 , B4 , B8 , then  add the  appropriate components  as per  the  binary  expansion  of the  exponent.  In our case, we add B8  × B2  × B to get the final product. It turns  out that these constituent powers of B (i.e.  8, 2 and 1) can be computed  in a recursive manner quite efficiently. Here is something  from the following link that computes  nk  for integers  n and k.

1:  int power( int n, int k )

2:  if k == 0 then

3:       return 1

4:  end if

5:  if k is odd then



2
 
6:       return (n * power (n × n,  (k−1) ))

7:  else



2
 
8:       return (power(  n × n,  k   ))

9:  end if

A similar approach can be used to compute  Ak   for a square  matrix  A.  I leave the nitty-gritty details for you to figure out.  Keep in mind that with NEWMAT on your side, the matrix  operations are structurally the same as that for scalars.

 

 

 

 

Complexity of Repeated Squaring vs.   Brute Force Multi- plication

 

If you have to compute  Bk , and assuming each multiplication is O(d3 ) (nothing clever here,  straightforward matrix  multiplication), and  since we have  log2 k-

log10 k

many  of these  to  do  (or,  since  log2 k  =


log10


2    we  can  replace  log2 k  with

O(log k)),  this  entire  operation   takes  O(d3  log k).   This  is in  contrast to  the

O(d3 k) procedure  for multiplying B  with  itself  k times.   Resulting  in a total complexity  of O(d3  log k).  If we used Strassen ’s method  for matrix  multiplica- tion we would have a procedure  that is O(d2.81 log k).

 

The Programming Assignment

 

1.  (Using NEWMAT ) I want you to write a recursion  routine

 

Matrix repeated squaring(Matrix B, int exponent, int no rows)

 

which takes  a (no rows × no rows) matrix  B and  computes  Bexponent   using the Repeated  Squaring  Algorithm.

 

2.  Your code should be able to take  as input  the size and exponent as input on the command  line.  That is, if we want to compute  Bk , where B is an (d × d) square matrix, I want to be able to read d and k on the command- line.  It should fill the entries  of the matrix  B with random  entries  in the interval  (−5, 5)1 .

 

3.  The output should indicate:  (1) The number  of rows/columns in B (that is read from the command  line), (2) The exponent k (that is read from the command  line), (3) The result and the amount of time it took to compute Bk  using repeated squaring,  and  (4) The  result  and  the  amount  of time it took to compute  Bk using brute  force multiplication. A sample output is shown in figure 1.

 

4.  I want you to provide  a plot of the computation time (in seconds) for the two methods  as a function  of the size of the matrix. That is, I am looking for something  along the  lines of figure 2.  For this,  you will have to place timer  objects  before and  after  appropriate portions  of your  code and  do the needful – as the following lines of code illustrate.

 

time before = clock();

B = repeated squaring(A, exponent, dimension);

time after = clock();

diff = ((float)  time after - (float) time before);

cout << ”It took ” << diff/CLOCKS PER SEC << ” seconds to complete” << endl;

 

1 See strassen.cpp for ideas.

 

 

 


In terms  of the  value  of the  exponent, tell me the  regions where one al- gorithm  performs  better than  the  other,  as  far  as  computation time  is concerned.

 

 

 

 

Figure  1: A sample output for different exponents  and matrix-dimensions.

 

 

3

 

 

Figure  2: A comparison  of the computation-time  (obtained experimentally) for brute-force  exponentiation and  the  method  of repeated squares  of a  random

5 × 5 matrix  as a function  of the exponent.