Programming Assignment 7: Repeated Squaring Algorithm Solution

Later on in the course, when we deal with the problem of computing the price of an European Option e ciently using the method of Dynamic Programming, we will have to deal with the problem of taking an n n matrix A and raising it to a large number. That is, we will have to compute Ak for large values of k.

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

    1: C=A;

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

    3: C = A*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(n3k) steps, as each matrix multiplication is O(n3) and there are k-many of them.
A candidate algorithm that is more e cient than the one shown above goes by the name of Repeated Squaring, and it described below.

Repeated Squaring Algorithm

Suppose you want to compute A11, 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 A. That is, you compute A; A2; A4; A8, then add the appropriate components as per the binary expansion of the exponent. In our case, we add A8 A2 A to get the nal product. It turns out that these constituent powers of A (i.e. 8, 2 and 1) can be computed in a recursive manner quite e ciently. 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
    6: return (n * power (n  n, (k 2 1) ))
    7: else
    8: return (power( n  n; k2 ))
    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 gure out. Keep in mind that with NEWMAT on your side, the matrix operations are structurally the same as that for scalars.


1





Complexity of Repeated Squaring vs. Brute Force Multi-plication

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


many of these to do (or, since log2 k =    we can replace log2 k with

O(log k)), this entire operation takes O(b3 log k). This is in contrast to the O(b3k) procedure for multiplying A with itself k times. Resulting in a total complexity of O(b3 ln k). If we used Strassen’s method for matrix multiplica-tion we would have a procedure that is O(b2:81 ln k).

The Programming Assignment

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

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


which takes a (no rows no rows) matrix A and computes Aexponent 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 Ak, where A is an (n n) square matrix, I want to be able to read n and k on the command-line. It should ll the entries of the matrix A with random entries in the interval ( 5; 5)1.

    3. The output should indicate: (1) The number of rows/columns in A (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 Ak using repeated squaring, and (4) The result and the amount of time it took to compute Ak using brute force multiplication. A sample output is shown in gure 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 gure 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();

di    = (( oat) time after - ( oat) time before);


cout << "It took " << di /CLOCKS PER SEC << " seconds to complete" << endl;


    • See strassen.cpp for ideas.




2





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 di erent exponents and matrix-dimensions.


x 10−4

Computation time vs. Exponent Plot

3.5






3






2.5
Brute Force Multiplication



2






1.5











Repeated Squares

1






0.5






0
50
100
150
200
250
300
0









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.








3